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The 2leorithm ASPEX is used for the 
dimensional plot with an azimuth, of 
alttaitade: of 8°) from the horvzomtal 
center of the graph to the observer 
drstanece, Ol azo. kmn from the -eraph te 
observer. 
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The main reflected and head waves for the 
model an fieure 20. Three generalized rays 
area tsitown O@isied vet eCP)ns aie, ava0: eC Pye: 


=2, Va0 (CS)). 


Synthetic. Setsmoocrans for a vertical: point 
force tand ‘Wwania tions ron ithe moded> shown) .in 
firure 20. The only parameters! that ane 
changed are the P and S velocities in the 
halé space.’ The P velocities: vary trom 5600 
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jae Introduction 


ine propagation of \elastierwaves ina layered solid, 
due to point sources, has been studied extensively both 
Gheoretically and numerically. ‘The’-solution. ether in 
exact On im approximate form, has been eiven for ditferent 
simple and complex models. My approach to the problem 
involves generalized ray theory and a Cagniard-Pekeris 
inversion because the method yields the complete and exact 
solution. 

The problem of an SH-torque situated inside a hayer 
overlying a half space was solved by Pekeris et al. (1963). 
iherasis “Of ithe ‘torque was vertical and the displacement 
dtc 0 Ehe source had aA -triangular shape at large distances: 
TMros sey pe O1 jSOurCce [excites only Horizontally polarized 
(SH) waves and these retain their shear character upon 
rerlection at the two boundaries. Abramovicd . C1970) pre- 
sented the solution of a compressional pulse in a layered 
babeouspace. =“ [hervertiical and horizontal components (om (tne 
displacement were given analytically in terms of individual 
generalized rays. im has paper, Nhe extended the work ‘of 
Pekeris et al. 101965) by considemine the ceneral case for 
Potssom eoratic © and nol the pamrriculsrecase(o7= 0.257 for 
the layer. He computed synthetic seismograms for various 
dephhes of “the source. and for different locations of the 


receiver. 
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Abramovici and Gal-Ezer (1978) studied the motion 
Otgasscolidgmade hupeofpsoneylayér over parhakiaspacecduejto 
Bhenprésencesnot aavebtiecal point forcejactinesberneath or 
attthée freessunfacesusing the generalized, raysexpansion: 
Each term in the expansion represented a group of rays 
beklectedrthexsameanumber of» timespas\PeoregS« «The time 
dependencesofisthesapplhied,forcerat thessource had: a tri- 
angular shape. Theyumodetled the vertical force as a 
SLKCSSAGaSContinudstytabong?they vertical Zeaxis and, they 
catoutatedtthenvertical and-horizontal displacements.for 
agreceiverriocatedsatethe freersutface. 

In this thesis synthetic seismograms for the general 
easerot Landhortzontallarb@irary i iorcesin’ therx-yfpleanés»are 
computed using the generalized ray theory. The theoreti- 
cal development of the problem was solved together with 
BrokescsorukedAbramoviciedurinel hisivistteon agsabbatical 
leave in 1979-80. The algorithms for the numerical 
evaluationnwerescarraed out: jointdyests they University of 
Alberta and at Tel-Aviv University. The x-component of 
the force gives rise to vertical and radial displacement 
waitevyéthety-componentegiv ¢semaseytesgtransverses mots on-. 
The complete mathematical solution for the case of an 
herizontalvarbitrarye pointet opcensatuatedsinside”™ a, layer 
overdyding acohadfoespece! isagivengin| Appendix, Ay which ds, an 
expanded version of the paper by Abramovici, Kanasewich 
and Kelamas? (498 20% 
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solution for particular boundary-value problems the solu- 
tioncplorsalli three sources, (SH-point tonuque, vertical 
Voane foOrece-and horizontal point. force). weresexamined. in 
detou Neen pnceiattenatiecal- solution is exact and is obtained 
inea elosed anailytical-f£orm, ashapsuperpostti onp_ofeindivi- 
dGualweeeneralizederays.waAllothéetpassibhekueysrwith surface 
and head waves are included in the seismogram up to a 
DigeLcularstimewof hint enes ts 

Another.innovationscot,this thesiseis thepintroduc-— 
ELOUSOL Jattenuation .andvdispenston im theesynehetic seis-— 
mograms obtained using the generalized ray theory. This 
<r Gene. iby decomposinge the selsmogram into individual 
sceneralized trays and incorporating the elLfects oct anelas— 
ticity in the frequency domain into each ray separately. 
The linear (theory of viscoelasticity alone with Futterman 's 
model is used to model the anélasticity and reflection 
coefficients for anelastic media are calculated to take 
LEO waccount ihe eLrect of the viscoelastic interrace:, 

im thas chapter ay review of the literature. about the 
attenuation of selsmic waves istmade. Particular emphasis 
is given to the frequency dependence of attenuation and to 
thesvarious mode Ls of -absorptiongeand, dispersion. In the 
second chapter synthetic seismograms for a thin "weathered" 
bayer) mndmeor a, Crustal \layer model ave presented anda dis— 
cussed. The third chapter includes synthetic seismograms 
Withethe wettecte lor attenuation sand dispersion. "Final ly-, 


inmvthestourth chapeerssynthetics are ‘computed for an 
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tion of wide angle seismic reflections and head waves. 


ieee etrerent derinittons of 0 

It is common experience that as a wave propagates 
through real material, wave amplitudes attenuate as a 
pesulthofsasevarhety, ofsprocesses,which+scan /be;summarized 
mMaenoscopically:as internal wirictiongor-anelasticity., The 
Deeustudepotieinternal,friction is expressed in; terms’ of 


the dimensionless parameter Q which is given by 
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where AE is the amount of enerey dissipated per cycle for 

a harmonic wave and E is the peak Strain energy in a cycle 

Of, tne harmonire Wave. “Physically; 0 1s a measure of the 

imoerreect lone 10 the edaSstLcity Ofbuthe. materials. 
Thesapplicapilbity cl thevabove (definition is sather 

limited; more commonly one observes either (i) the temporal 

decay of amplitude in a standing wave at a fixed wave- 

number or (11) the Spatial détay*an a propagating wave at 

a fixed’ trequency. "AsSstming thar @ 22>) and for’ media with 

linear stress-strain relation we get that the wave amplitude 
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teis tnesangular frequency. Thexsquantitysi/Osis.also 
freterreda tc as the loparithmic decrement 6. YeEduatidon (1.3) 
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Pres Jscdetiped sing Eq. i144) te called, temporal 0. 

im iecase (il) |one Lodlows,) a particular wave peak alone 
a distance dx and.observes the» gradual spatial decay of 
amplitude. Assuming that, the direction oft propagation ds 


a1S0 then cirecttonyol, attenuation wesycer 


Wasa (21255 
dx 
Wheres is the wavelencth. JUSine stbhne elation petween the 


wavelength A, the angular frequency w and the phase 


velocity c 
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Wlon thetobvioussesolution 
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A(x) = Ale E9 (2.8) 
where AY isethepneampligude at x = 0.” Equatdon (128) defines 


Ehewvalue of ‘the spatial Q. 


The absorption coefficient igevsis—defined as 
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Equation (1.10) describes the attenuation model that 
wilintbestsedwmin thisethesis, Paeaggannexponential decay ros 
amplitude with distance. The above definition of the 
absonupttongecoersicienteimphiesmehatucertsyarhanearyftunction 
of gdréquency¢eaemore detailed discussion of this topic will 
Occup in later sections: of this chapter. 
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The ,éLiniéeiontoltOpeiven@in ediarron*€1l.1) is 
equiva bent eto thhe 400i Can oseilLlavorytelectricalVerretit. 
TrRiseprovides «an talternative ‘definition *for*Q which référs 
to the forced vibration of a specimen with resonance 
frequency. ni. Litthe tfredden¢ytofithe vibration is varied and 
Les anplttude. held Constant, then, the tanplitude ‘ol vibranion 
of Lhe specimen traces vOut «a2frescmanace curve’. By measuring 
the width Af of the resonance curve when the amplitude of 
vibration equals 1/2 of the amplitude at the peak frequency 


f, Q is expressed as 


Af 
ue he eect 
F Chl. 


O|r 


All these different definitions for Q may be applied when 


GO) =e)il Vand they imply that the direction of propagation is 
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inMe=sane-as* the dirTrettion OffSttenuati on (Zener.9 1948) . For 
strongly attenuating media, these definitions will disagree 
and Sthus*werhave to’ detinerOuditicerentivys @ 0" Conneda® and 
Budiansky (1978) proposed the following definition for Q, 
which is probably the best of several different definitions 


that scientists have used for Q. 


Ae een thy 

0 Re[M] CAG.) 
where, Die hg ae eZ for P waves 

M= u for S waves 


A, HW are the complex Lamé parameters. 
Boer apove derination for 0 will be. wsed for calculating 
fer lection and “transmission ‘coeriicients for anelastic media. 
In summary, the common measures of attenuation are 


Genrineda sand telacecd. to one another as toliows: 
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1-3 Experimental measurements of attenuation 
The measurement of attenuation for earth materials 
is well covered in’ the geophysical literature. The values 
of Ocanrbe dervermined in) the Laporatory,. 2nvene freidyor 
by analyzing earthquake and normal mode data. In this section 
a brief review of the experimental measurements of attenuation 


will beemade considering only Laboratory data and Eield 
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experiments. Attenuation measurements from earthquake 
records and normal mode data will not be discussed although 
Bie@hite patune fis ,.quite extensive inethisvaread. »kRarticular 
attention will be given to the frequency dependence of Q 
and) the tlineaygityrokuthetabsonpthion coekiictient rin terms of 
HReGUucn Gy ayeBefore yooinguintostheenaingdiseussion,one has. to 
realize that the accurate eaeunenenn ofLsattenuation wis a 
difficult task with experimental and interpretational prob- 
lems. The choice of a particular technique is based largely 
on the frequency range of interest, the actual values of 
attenuation jpandythesphysical»,conditions sunder which the 
Samplevwill bexstudied. The methods generally used for 
measuring attenuation in the laboratory can be classified 
intonehe stolkkowing categories +( Zener, LO48esKoleky, 19:53 
penbpetbear set dalajiL97 3): 

he Free vibration: Inathis,techniquesa.~rod of srock 
is «suspended vertically so that a) mass with.«a, large moment 
Ofpinertiagcanabe rattached ttoutts glower pends «hiathe mass uis 
enventakkick:theesystemeyibsatespireelyce EBhe ratenol decay 


of«the amplitude of \these, free,oscillations jis ,attributable 


ie enercyehoss Lngthesrock sp rDefine theglogarithmedeerement 
as 
gn(A,/A,) 
Lth2 
One ea ey Cie 3) 
(t,-t,)f 


where Ay and As are the amplitudes at times ti and to ad af 


is the natural free-vibration frequency of the system, the 
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Vorwver Or Oils calculated Usineged /Ov=" 677s 

Zs Foreed vibration: ‘The’ method’ has been mentioned 
already where Q may be found from the width of the resonance 
eurve. Thquetron Cl to) WVerased in thes calculation of 0% 

3. Wave propagation method: The transmission of a 
Pemsewtucouch the rock Vand its" detectron te stthe basis of 
this technique. The method assumes an exponential decay of 
amplitude of the seismic pulse with dvstance om tinve=and that 
ere wcans correct Lor Losses ‘due ‘to’ spreading’. © etiections’, 


veer ace Tons, etc. 


ie VOUSerVaAtLOn Oh uStressSstrain curvecs = tor cyoL peal= 


Peeoccescsed. Materials the enersy Loss sper “cycle,” AES" i’sh'the 
area of the mechanical hysteresis loop, in which the phase 
On the strain lags benind that of@the Stress.” “Thus, O° can 
pewealeulatved: using equat Lom “Cl. 19% 

The use of stress-strain curves will provide informa- 
Clon at frequencies below 1) Hz, white resonance’ vibrations 
are used to measure the attenuation in the range of 100 Hz 
to 100 kHz. Wave-propagation experiments are commonly 
restricted tte the Uwltrasonic range ot 00 Kg “or hieher.” «Let 
us now discuss some key papers dealing with laboratory 
measurements of attenuation. 

Birch and Bancroft (C1938) measured ‘the’ 'O of ‘a ‘core 
Of Quincy granite using the resonance technique They 
concluded that Q of Quincy granite was approximately indepen- 
dent of frequency in the range of 140 to 1600 Hz. 


Born (1941) used the same technique to measure the 
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orarithmic decnement, ot iseveralirock joores takem from wells. 
Pes ound that the logarithmic decrement) forthe drwy rock 
samples was independent of frequency. His studies and 
measurements with water saturated samples showed that the 
focal decrement was the .cum)iofs two terms... the, first, term was 
frequency independent and the second one which was 
Peoport tonal; to frequency. 

Bruckshaw and Mananta (1954) measured the quality 
Pac tor won Six ditierent rock cores Ceranitte, “dolerite, 
diorite, sandstone, shelly limestone and oolitic limestone) 
Pen c si requency, range 40, to, 120 Ha. They sconcluded that, .0 
was pEactically<independent, of ~trequency, for these  wix rocks... 

Peselnick and Zietz (1959) used a pulse-echo technique 
to measure the absorption, coefficient..in three different 
limestone samples. Using a-~piezoelectric crystal attached 
Eo one vend of the specimensa short, pulee of sonpressional or 
Shear waves is, cenerated.. The first arrivals and, multiply 
rerlected yuenels are detected at! the tother jad Got the 
Specimen py a second Crystelie) Lnese (arrivals are displayed 
on an oscilloscope and) the decay, in amplitude, is, measured, 
thus me tne. apcorpt tony COCt ia Client scan se Cael ated. Js) thew 
found tnat the, absorption jcetieileng was) proportional, £0 
thes prequency Gn ithe range of 3S hateo elo Miz for compressional 
and shear waves. Tieabsor pion icoenni cient caleulated this 
way ineludes, the, effects of apparent, losses due to geometri— 
cal spreading and the reflections and transmissions at the 


possible specimen boundaries. However, they found that these 
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bosses were meglicible an’ respect to the: absorption. 

Wye teal CLI62) stud ted. the apsonptLtonnoLese.amic 
waves in fluid saturated porous rocks. They proposed that 
ene logarithmic decrement was the sum of two terms. The 
first term was practically independent of frequency while 
the second term was frequency dependent. They studied the 
abSOrption in terms of ;sfluid saturation and they did not 
give any evidence showing the frequency dependence of the 
total logarithmic decrement. 

Knopoff (1964) summarized the laboratory measurements 
Ope Jin metals, non-metals and rocks. [pelea neLowe wwe tn 
Ene schnortest title in the. geophysical Literature, ~Q", 
Knopoff indicates that for small strains the absorption 
Pocwricrent as 4 linear Tunetion of frequency sand hence © is 
constant. 

Bradley and Fort’ (19660). Atwell and Ramana (19.616,) 
reviewed the literature and they concluded that Q is inde- 
pendent of frequency over a range of about 10° HZ Dey 
tabulated laboratory and field measurements and measurements 
from earthquake records. Atwell and Ramana applied a least- 
SquUaLes analysis wand concluded Chat Ene absorpbion coeiitietent 
Was a linear function of frequency 1m the range Hoe Hz to 
LOG Mirz. 

Toksoz et al (1979) studied the attenuation of seismic 
waves ain dry and, saturated rocks. Their laboratory measure— 
ments were made at ultrasonic frequencies (0.1 to 1.0 MHz) 


and they showed that attenuation coefficients increase 
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linearly with frequency (constant Q) for both P and S waves 
dnebothedrysand fluid saturated rocks. They found that the 
attenuation in water saturated rocks is greater than in dry 
rocks;+while»the, attenuation’ in, frozensrocks is very»much 
lower than in saturated rocks. They also studied the effect 
of pressure on the attenuation of seismic waves and they 
concluded that attenuation decreases (Q increases) with 
inereasineypressurestoreboth,PtandsSiwaves inealilicasesxof 
Saturation.y Phearaterofpinerease,of .Q se high at low 
pressures and levels off at higher pressures. 

Therefore,g-in lLookingijat the daboratoryydata,which 
besa justebecnadisccussed awewtind thatpiortgdryeanocks andtmetals 
the absorption coefficient depends on the first power of 
frequency and thus Q is frequency independent. Laboratory 
measurements with fluid saturated rocks do not lead to the 
same conclusion. iIneltquids a istproporGional tomitrequency 
(Pinkerton; 1947), so that in, some highly porous: and, permeable 
rocks the total ora may contain a frequency dependent 
eomponentP(Bernes! 941s Wylie et ‘al, 1962). ° There is also 
evidence that on is independent of frequency for saturated 
rocks as pointed out by McLeroy and DeLoach (1967) and Toksoz 
et al (1979). Most authors suggest some form of frequency 
dependence) for 0, foryvtluid saturatedprocks.g#This formananges 
from Oma vatying,dineetly asy)irequencysto pee varying 
inversely de the square root ‘of frequency. Ft ic evident 
that the rock Lexture ,-ponosi£ty, pressure, Cemperature and 


type of saturated fluid play an important role in the 
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relationship between attenuation and frequency for fluid 
Sacubated rocks’. 

Let us now discuss briefly the field experiments 
nu -situ) tor measuring 0. If the objective is to under- 
stand the attenuation of seismic waves in the earth, then 
Eremamportance Of measuring/0 in situ is) apparent. The most 
common method of measuring the attenuation in situ is using 
Carawepom vertical seismic profiles .VSPjIe Seismic pulses 
generated by a surface seismic source bey recorded by a 
single downhole seismometer positioned at various depths in 
a well. The lateral separation between the source and the 
well is generally small compared to the seismometer depth. 
The spectral ratio method is then applied and Q estimates 
are obtained. The spectral ratio method involves a Fourier 
atvaryeis or the date from which we"can compute the value of 
Qsand tire phase velocity. “lo explain thesmethod consider 
a pulse which has travelled a distance rT from the source 


and another pulse which travelled a distance =) iron tive 


where Fi fw), F,(w) are the complex spectra of the two pulses 


and k is the wavenumber. Make the wavenumber complex 
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where c(w) is the phase velocity and o« is the absorption 
eoeftfiicient. Substitubesine(lseLo) landerakevVtheretntwoteboth 


sides we obtain 


-a(r -r,) (alee? a) 


and 


w(r,-r,) . 
e(w) = - $5 (a) -4, Co) Cis Zy) 


where A_,A 


Ao are the amplitude spectra and b1 25 are the 


Pics e spec ura, of the two pulses.” Equations, Gliu2.) and (Cl. 22) 
define the spectral ratio method. EGudelon -Glecd), indteates 
Pate pO, Ol rive ina tural: logarithm jor the ratio, of. the 
amplitude spectra against frequency can be used to determine 
the nature of the frequency dependence of the absorption 
eceriicient. Equation €1.22) defines the vdispersion, curve., 
woe. tieophnase Velocity asiva functione ore rrequency. 

The first such experiments were carried out in the 
Prerre Shale area near Lamon, Colorado.” This particular 
location was chosen because it provided 4a thick section 
CAapproximetrely L000 m) of unitorm shales veny close to the 
Cartchies. cur face. 

Ricker (1953) performed a series of field experiments 
inthe Lierre shale area in order mo test the validity of 
his wavelet theory. According to Ricker a modified wave 


equation which includes a term with a first derivative with 
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respect to time, is used to explain the seismic phenomena. 
Ricker's model leads to an absorption coefficient which 
Varles@as®thessquare ofvtherfrequency andSit predicts that 
taecewavelet, Dreadthrwilivincreasetashthe Squares tooteoflthe 
exaved™time. SHe presénted' graphs  whith® show the validity 
becuse “neoryt “However, daSthe®discussion of his paper many 
authors pointed out serious problems in the comparison of 
observed and theoretical wavelets. 

McDonald et al (1958) made attenuation measurements 
in the Pierre Shale area, in the same area where Ricker did 
HLSUworks Shots. at-depths" of 9764 2°to89l . 4Gemiwere,fired®into 
Pives peophones“at> depths*from’ 1O79to 229 me ~They®used 
velocity log data to apply corrections for geometric spread- 
ine and they’ measured the’ absorption by Fotrier analysis of 
the corrected waveforms. They: coméluded@thaeswtioreverticatly 
travelling P waves the absorption coefficient was a linear 


Punct1 on oOfetrequencys® invtnrerrange of 100 to 000) Hz. . They 


obtained@samilare results. Por’? horizontal travelilanges waves. 


They also studied the decay of wavelet amplitude with travel 
time and they found that the wavelet broadening with travel 
time *was Less*than Ricker had measured 

Wuenschel (1965) studied body wave attenuation and 
dispersion. The experimental data used in his study was 
the®pulse=propagati on® experiment? donee byh Mobidtdine the’ Pierre 
Shalet X(MNeDonald? ct ai) 1958)° and) as model seismic analog of 
thee Mobidlexperiment wsing ultrasonic’ pulses in a® plexiglass 


sheets “Her showed that the’ dispersion phenomenon was’ present 
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along with the observed absorption and he concluded that 
thes ¢2edd. experiment, conducted: by. MecDonald,.cet. alp.(1958). in 
mice Penne Shalejicondinms Fatterman's theory €1962) . 
Futterman's theory provides the attenuation model used in 
this thesis, and a more detailed discussion will be found 
bngsohe nextyjsecitdons* of: thissichapter. 

Cole (1965) studied reflections from the Newfoundland 
abyssai plane and.the Gulf of Alaska. He approximated the 
reflection process. by a) three-layer fluid model, with 
aeeenuarion sin the second and. thind, layers... » He, then. calcoul-— 
GuedTheonmetrcal reflection, coefficients considering the 
Segiment- attenuation. as. proportional to. the one-half, first 
and second powers of the frequency for comparison with 
measured coefficients. For the frequency range 100-1000 Hz 
hee founds that. the; firstrpower frequency attenuation, Laws in 
the sediment is consistent with the observed dependence of 
the ocean-bottom reflectivity. 

Tullos and Reid (1969) measured the attenuation of 
Seismic. enereye in- the» sediments. inp the Gul, of Mexico... They 
used. .byias tdng, icapsi as sources; and, peophenes: cemented sto) the 
Aarth ats vardous. depths dna borehole, as uneacedwer she a they 
applied\correctionsr to the. data due-to geometirical spreading 
and theys memoved) the effects, offpconstructive and destructive 
interference. due. to: reflections by. averaging the, spectra of 
hany. traces; over a.small, interval in depth... The, ratio. of 
spectra at different depths was then calculated and from 


this spectral ratio they determined the nature of the 
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beequency dependence*of the absorption: “They found that the 
absorption coefficient tis atlinear "funetion-of "frequency “in 
Biveirange .fromes0 te 400 th2. 

Hamilton (1972) studied the compressional velocity 
aud Yabsorption*in the “sida (floor “off San Diego%in ‘the frequency 
macaw. wo mkiz eto VLO0 SKHZ Me tconelided Sihat vattenuation is 
approximately dependent.“on "the first power "of frequency and 
Siaceve PocDiy dispersion is almost absent in water saturated 
sediments. atte also indicated “that trtercrain tirtct ron 
appears ‘to be), “by far, the dominant cause “of wave energy 
damping in marine sediments. 

Ganley and Kanasewich (1980) measured the absorption 
and dispersion from ‘seismic check-shot data using the spectral 
ratio method. The data was taken from a sedimentary basin 
in the Beaufort Sea. They showed that frequency dependent 
itesisesi wie tho reflections ‘and *‘transnisstons play “an “important 
mo Nelitands tiie yetapp 1 edi tcorreectonse dive vo these “osses using 
synthetic seismograms. Their analysis showed that the value 
Om tOlkis Talmoct® frequency independent and that the dispersion 
Messina any bie dat a? ais! Consi's terme: wreuhn Pubtternan s model”. 
Measured? Qf values™were! 4322) for, a depti> tnterval from 549 
Hoe IOem land 6G726MLor a depth anverval arom 94 )>-*vo Slim. 

Hawoe: ( 1osdpetised, data’ Eromvertical™s eirsuic profites 
tol meneunre: Vetenuation vie Mie spectral” ratio method." The 
daiter werel taken? ‘trom five detailed velocity surveys, one was 


carried out in West Texas and the other four in the Gulf 
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Coast area. Actual measured values of attenuation varied 
Dy a factor or !0 "ranging ei rombless than tOUDCton0 sod B oper 
wavelength, depending on the lithology. He found that porous 
sands have much larger attenuation than neighbouring shales. 
studying” the shape of the seismic pulse, he concluded that 
attenuation changes the pulse shape noticeably, which is 
ancther™“indtecatton that,.dispersion is dan essential consequence 
SLPapSsorptions 

Spencer et’ al (1982) studied the basic’ problems 
encountered in extracting estimates of seismic attenuation 
irom data recorded*en3verticalvseismiteprofidiess tiheyeaepaid 
Pacticular attention to interferences etieets, bepatial 
resolution and frequency dependence of attenuation. They 
found that’ for -smalil*receiver’ separationes thevattenuation 
computed from the spectral ratio method is much more strongly 
intiguenced by the*® Vocal ‘stratigraphy ine enerdinumediateswa cind ty 
of the seismometer than by the attenuation in the depth 
interval between seismometers. They also showed that the 
spatial resolution is strongly! influenced? by) thei Local 
stratigraphy and°in most’cases attentuation @stimates sroudidunot 
be possible to'awdentiiy Litholopivtes: Or tcond.t loniss which are 
the source of anomalous dissipation\ Pa ott ines ther amplitude 
ratios @eainst frequency. they Hound Gharevehe teraphitexhibits 
a linear trend which is in good agreement with previous 
experimental results. 

Summarizing the field measurements of attenuation it 


is evident that the absorption coefficient is a linear function 
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Stitrequencysend) hence Os isi Const ante This statement is 
Emuerior frequencies: ing they ranges of L0.toiS00,HzZ ands for 
Pyecituycocksa)yVerticalsseismicyprofides! (SP) eprovide the 
most successful method of obtaining in situ attenuation 


estimates. 


1.4 Attenuation mechanisms 

In order to evaluate and interpret laboratory and 
field measurements of attenuation, possible attenuation 
mechanisms involved are needed. Numerous mechanisms have 
been proposed and each may be considered to have a greater 
degree of importance to the overall attenuation under 
Cextainy physical conditions. in this Gectiony ay brier review 
of possible attenuation mechanisms will be made, considering 
only these mechanisms responsible for seismic attenuation 
inlupper crustalerocks) Studying attenuation mechanisms one 
has to keep in mind that each mechanism depends on rock type, 
SaAkurationystates, pressure,.frequencygranges amplitude of 
thesacousticvwaver andwotherjwardous tock propertiess 

Walshe({i966) proposed, thatethesmain>sourcesofaattenu- 
Ationtingdtyt rockswis+basedsonathenitictionalijdissipatdongas 
Cxackpeureaces on contact slide Yrebative’ to one another 
during passagerofaaucseismicgwavesy, Hegeonsideredgthateald 
other sources of attenuation are lumped together in the 
intrainsael attenuation, »which;mighttinclude lossesrdteyto 
thermoelastic effeckssiviscosityseete. inphiseanalysisse 


Walshyv formulated’ the» problemuby approximating) thex»crackstas 
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ellipsoids and then he calculated Q values for P’°and S waves. 
The resulting expressions are rather complicated, neverthe- 
less, the Q values obtained are frequency independent. 
Johnston et al (1979) utilized Walsh's formulation to model 
pressure “dependence for “the “attenuation of sultrasoniic waves. 
They showed that increasing pressure decreases the number of 
efacks Contribpubing (to as etenuationebhy@inricbion,+ohus the 
attenuation is decreased. Other authors (Mavko, 169. 7-9 8 
Miter zetval, vb979 \°proposedethatestrietiontistnot ea dominant 
mechanism particularly at low strain amplitudes and low 
Frequencies. However, under certain laboratory conditions, 
namely high amplitude ultrasonic experiments, it may be 
dominant. 

Savage (1966) applied Zener's theory of thermoelastic 
attenuation to explain the attenuation of elastic waves in 
solids. He showed that thermoelastic losses are associated 
With toneeftlatltecavitiestwhich* @an réeptesent cracks sintehe 
medtunty Thistmode? predicts tatdecréasexin att éenuationswaith 
inéreasing i pressure and ancinctreasestin @fforelownhinequencies. 
Recently, Armstrong (1980):\"Has*proposed another thermo- 
Gelastic move! foe which attenuation “is essentially frequency 
independent. 

Ina series of papers, Biot 1956, 1962) developed a 
mathematical theory for the dynamic response of a linear 
porous solid containing compressible fluid. According to 
Bi6e' setheoryeartentation*ig aorésudltaot sthevmotion!ofethe 


pore fluid relative to the rock frame and it depends on bulk 
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pOck properties: “This typeret mechanism is generally con— 
eluded to produce negligible attenuation at) low frequencies 
pueeonsolidated rocks "(Whitey hl965).% Howeverg™it maybe 
important @at = ultrasonic#irequencies™(JOhnston et alkkov979) 
ore inv’ permeable; unconsolidated *sediments at sintermédiate 
frequencies. 

BfOt-type mechanisms considerstoniyvbulkfiluidt flow in 
Porous Fock and they ignore inter= and intra-crack, flow 
(Meocalerlow) {Mbothtot4whith maytdissipate sed smi cwenersy < 
Pmeercrack’ fluids filows “sometimes®known. as Usquiretetiowylwas 
first proposed as an attenuation mechanism by Mavko and Nur 
C1995) % 

Another type of mechanism is a stress induced 
diffusion model based on the thermally activated motions of 
aeonowor defects on the lattice of a crystal ‘under the 
ime lueneelonm-externalvetresses (Tittmann etal, 1£980):. This 
type, OfhkRechantiesn is Similars rer thatnduer bot vistous* erain 
beundary\dampine, @wheré! there Gs*airelaxation of StYressevat 
grain boundaries. Both mechanisms are relaxation mechanisms 


and they are generally described by 


! 
21D = oe 
ltw t' 


O|H 


wheres®D» ais’ the ratio®or the? createst® non=-elastic stress to 
the edasctic Srvain,*t! @s atrelaxation time*which varies 


with® cemperature Taaccording®*tofan= equation of*the form 
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where G is an activation energy and k is Boltzman's constant. 
Summarizing the attenuation mechanisms it seems that 
the *most important absorption mechanism for sedimentary rocks 
@ieseiemicyirequenciestwould be*slidinge.at cracks or grain 
boundarpes Se 0This emechanisi “impli és "atphtase “anele between 
Serres tand Sirain which Leads “to va hysteresis loop. Relaxa- 
Pron tnecwayisme such Vas “rain “boundary «damping * “thermal 
CUECENGS “Or stress induveéd difiusion play -a*secondary “role 
in the attenuation of seismic waves. Finally, the presence 
Wawa Ehibutd evs: G4a*third *possible source "ofr “absorptrfon?’ ” "Brot— 


! 


@ pe techianmens *and Visq@ire Sf Low "can "be Mconsidered tas thie 


most \donimant Mechanism sin “porous “rocks: 


1.5 Attenuation models 

Dine tie cece ty) “f or. Ginidorporating sthe ve Pirecu of “att en— 
UtaALoOm Onppwrane. "propde ah on siniYrealrstie “media “rs “obvious to 
every geophysicist. Several models and mathematical theories 
have) ibéen developed which take into account the effects of 
anelasticity. In this section some key papers dealing with 
attenuation models will wie *ditscussed. 

firettsimple st model “whieh "is "tire ‘basis low Purther 
generalization is based upon linear elasticity and Stokes' 
viscouniye 8th “are ritaie tot “thd s icind ts icalted viscoelastic. 
Riteker Ml 9537 aia cha s tcliasctica ll! paper tried ‘to model tire 


absorption by adding a single term to wave equation following 
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the work of Stokes (1845). He concluded that the absorption 
coefficient varies as the square of the frequency and thus 

Q is frequency dependent. Pes now Clear epbat Reekenus 
model is not- adequate to describe the anelasticity in earth 
materials since the frequency dependence of Q contradicts 
practically all experimental observations. 

Koiskyr (1956) ‘studied the propagatton of Short 
mechanical pulses along rods of three polymers, polythene, 
polysterene and polymethylmethacrylate. He first calculated 
theoretical curves showing the phase velocity and attenua- 
bEon £Or three models, the Maxwell model, the Voiet model 
and the standard linear solid. The Maxwell model consists 
Gtea pertectiy elastic spring in séries with a dashpot whienr 
fas Newtonian vViscosaty,,in thesVYoiet model the spring is 
felted sacross the dashpot.and.the standard linear solid is 
represented by a second spring in series with a Voigt model. 
He concluded thatthe standard lbanear solid model gives 
better approximataon to the behavioqur et fayrealivisco- 
Stachic Sseltd but sett) canvonly bermutted over, a timated 
frequency range. 

Futterman (1962) showed that dispersion is a necessary 
consequence of the medium absorption and is determined 
Ulambieuously by it. He assumed that the absorption 
Gociidecremt sc Ge) is striotly elinear in the frequency, over 
the range of measurement and that the principle of super-— 
pcsition, Ls valid. Consider a one-dimensional plane wave 


in the frequency domain as follows: 
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uCx,w)”>= Goethe C25) 


where K, the propagation) constant, can. be expressed in 
terms of the phase coefficient k(w) and the absorption 


eoefficient a(w) as follows: 


K(w) = k(w) + ia(w). Cit? G) 
Lhes., 
u(x,w) = homie Gua. eS (ates) 
oe x 
or wx) = Ne ee aS ety, Cie28) 
wiere  c = ¢c(w) vs sthe phase velocity. 
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Equation 1-29) Shows thatoifiithetabsorption) cocfificaent 
and phase velocity are known as functions of frequency then 
the effects of absorption and dispersion can be easily 
introduced in the frequency domain. 

Futterman assumed a low cut-off frequency wo? 
characteristic of the material, below which no dispersion 


exists, and he expressed the complex wave number KC) on 
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terms of the index of refraction n defined as 


TCO). = aS = Re[n(w)] + ilm[n(w) ] Ce 363) 
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where K 64) defines the nondispersive behaviour of K at the 
same “frequency. He then used the principle of ‘causality 
to derive relations of the Kramer-Kroénig type relating the 
dispersive part of the index of refraction of the medium 
BOet ne an sorplLive part by means of Tan intesral over the 
entire frequency range. These K-K relations are a 
comseqcuence Of the principle of causality and “follow without 
recourse to a specific wave equation. He considered three 
Gi erent LOrms 10% abSOrpthLon Which Satisty the “assumptron 
Of linearity with “respect "to Erequency. in taiws cnests 
Ere LAS ubPtenrman Ss 2abSsorptaon term te considered, namely: 
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Ceeniac1 on model used in this thesis. the Low, cutotLrte 
frequency Wo dOesi NOE COonrespond to-any -specaiic pwhysical 
mechanismrand it cans be chosen sutiiciently low. The walues 
CeO e sob tained "py ieguatton, Cl. scolreshowltnat 10 ic 
practically constant and frequency independent. Futter- 
man's model is in excellent agreement with experimental 
data (Wuenschel, 1965) and it has a major advantage in 
biimis We meed not. appeal. tothe physical details that’ are 
CUWenmacrenistic¢. lof a (particular theory. They are bypassed 
eo vreld the dispersion darectly from the absorption 
coefficient, in this case taken from experiment. 

Lomnaitz. (1957) trivede £0. exp laam the relationship 
between Cransitent: creep and interval friction in solid 
materials. He used, Boltzmann's superposition prineiple, to 


relate, stress and Stratn -as follows: 
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where y(t) is -theisrelaxation functions “Equations (1434), 
(1.35) pare-convolution integrals, thus. by: Fourier transform 
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Equation (1.41) is very similar to the frequency 

dependence cof (Q}givenSby)Futtérman inlequatione(l:33). 

Indeed, the two theories by Futterman and Lomnitz do not 


CeereLoecteniticantily ‘(Savace and O'Neill. 1975). Futterman 


assumes a complex propagation constant while Lomnitz 


Seonsidaenrs “a ‘complex elastic modulus, ~Neither theory provides 


any information about the physical mechanisms of absorption. 
Carpenter (1966) calculated the impulse response of 
a system satisfying the constant Q hypothesis. He compared 
the models of Kolsky and Futterman and he concluded that 
they are very similar. He also concluded that in many 
Gaseéseittisesmore convenientito*work inwthe frequency “domain 
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where c. is the phase velocity as the frequency becomes 
infinite. 
ptriaek CL967> 1970) ised the followine telation sor 


mremwapsorplLion coefficient : 


a(o) = ba (1.44) 


wieLe "Db ols a “constant and 0 <“s < 1. FRe then derived a 
dispersion relation for the above model and he noticed 
reasonable agreement with experimental data. He also 
developed an asymptotic expansion for the time domain shape 
Crea ibraveriline wave at Jarce distances: ‘He vused ja “solid 
of 0 = 50 and he compared his results with Futterman's 
theory. He found that the shapes of the pulses are similar 
but theytime ts very different. “The difference ain times 
Perdue tO svarlous approximations, that Strick made, “thas chiis 
theory 1s woalid only for, frequencies. much: greater ‘than, the 
low frequency cutoff. 

Kiartensson €1979) presented a linear model for 
attenuation with © exactly jindependent of frequency. 
According to his model the wave propagation is completely 
specitted by two parameters; e.¢. 0 and as a phase velo- 
€Ley Vat on aArbatrary reference trequency ee He started 
using Boltzmann's superposition principle with the following 
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He then derived a simple exact expression for the phase 


velocity 


as "a, Lunetion of; trequency : 
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K7artansson s model is causal since both the lereep and 


Tretaxsatiom Cunctrons Vanish hor) negative tame, no strain 


Can precede “applied Stress, nor Can anyestiress precede 


applied strain. 


Finally, the standard Linear solid model wiki be 


examined. 


This. model is well described by Zener (1943). 
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a -Step change in strain; P, is the relaxation time of 
strain under an applied step in stress and m is a constant. 
The Q is obtained as a function of w and the relaxation 
times 
GU) — 
oe cpg se) 
ane Sears , (1549) 
1+w e 
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The above equation shows that Q LS Obs COS tance DIbL et 
behaves as a resonance curve as a function of w. The 


aLeenuatron, O° ,. is concentrated near the :redquency 


ak . 
With O t behaving like w for frequencies below 


Vea 
this central peak and like et fOr Lrequenciles@apovie: 1 t.. 
The corresponding expression for the phase velocity shows 
AVMONGOEOn LCA increase With’ frequency. 

In order to reproduce the effectively constant ab 
Values observed at Seismic frequencies, Liu et al (1976) 
assumed, that attenuation is due to a superposition of 
different relaxation phenomena, each of which corresponds 
EO a relaxation peak and is represented by the Stress= 
Strain relation of equation (l.43) >. Thus, by superposition 
Of twelve telaxation peaks of thus Eype: they optained ia 
model with constant one Over Eneurance 00000 “to, LU he. 
Mey also showed that the dispersion relation for such a 


model is very similar to that arising from Futterman's 


theory. 
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Zedmiaintroducti.on 

Generalized ray theory was introduced by Van der Pol 
and Bremmer (1937) in a classical paper on electromagnetic 
Badiationraron aspoints source Thee lptesrads terms, indehe 
infinite expansion were described clearly as unique 
Pemtected or multiply reflected rays from the -boundaries in 
the layered medium. Beginning with Cagniard (1939), and 
peter ci( 1940). varioussyforms, of; this)technigque have, been 
appliaed. to the.study.of elastic, wave, propagation., 4 A 
detinitive.exposition was. given, by. Pekeris.et, al.(1965) for 
invel case ot a*solid layer over a half space. They describe 
the elastic wave solution as a double power series over 
imadees Tf and (iy o= OL...) in whieh -each ‘term, an) the 
expansion represents a ray that has been reflected uU times 
as Asecompressional.(P) «wave and 1) times) asy ae shear, (S) wave. 
The term “generalized reflection and transmission coeffi- 
cients" was introduced by Spencer (1960) and "generalized 
raver) by, Casteroas ety ad (1973.6 

In this thesis we follow the development of generalized 
ravi theoryatorethée. case) oft ay pointystxresss pulsexinyan 
elastdcelayer/over ay,halfyspaces; ,Exactisynthetic,seisno— 
grams are obtained by computer modelling for.a, number. of 
cases of interest in exploration and earthquake geophysics. 
These include the weathered layer problem where recordings 
closes fo the, sourcesmayy beg daft zsicultstorpinterpret.,sinee, they 
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field observations or plane wave modelling. Another case 
examines the simple model of a "crustal" layer over a 
"mantle" half space for a horizontal point force. Finally 
in the fourth chapter we examine the region where head waves 
Sepanave from the reflected arrivals in a typical sedimen-— 
fAarvyeceCe Lon consisting of a low velocity Laver jover a thick 
high velocity one. 

Consider the propagation of waves from a buried stress 
guscontinugity to a.isurface receiver in as solid ‘elastic 
Paver sOveron hali=-cpace., The exact: analyticesolution in 


the Laplace transform domain has the following form: 
CNG ERROR, Foran Dp Gyphbet Daye chalets ee C=) 


The transformed displacements (uy =. qs et V; re = w) are 
iMiecywlindrecalweoordinates: (rid..2)) An whic ci reutar 
symmetry is assumed about the z axis and the source is 
placed at a idepth d below the origin. The source function 

F will be a stress pulbsé Of type | witha trianeulaw shaped 
time dependence when viewed in the far-field. The algorithms 
wild be lcomputed for a horizontally polarized torque, (j=); 
a Vertically (j=?) and a horizontally, (i=3) concentrated 
force... The Laplace transform variable is p{=i0) where Ji 

is the angular frequency. M equals infinity unless compres— 
Sional waves ane absent in which case at is sero. The siear 
modulii, shear and compressional velocities in the layer 
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each generalized ray is) given by ey an integral involving 
integers wu and Vv, where yw is the number of times the energy 
is refiected as a P wave and v the number of times as an 8 
wave. the response’ of each generalized ray has the follow-— 
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where x is a scaled integration variable Gz > jin Witches 


is a spatial wave number measured in the direction of wave 
propagation. Je is the Bessel function, of order m=0 or 1 

’ : , : : ie p F 
and J Ls-itss derivative: R. and r. are generalized 
Pecvectien Ccoerticients for the 2th displacement for “the 
POtality Ol rays starting as P and arriving at the receiver 
after being reflected »w times as P and V times as S waves. 

: S) S Aes ‘ 
Functions Ro and r. play ra Sinirar *co'be "for crayss sit moans 
aero, AS Will-be ililtstrated, "each Pov pdr represents 
Several distince ray paths “and the generalized reflection 
eocetticients are, the sum of all their reflection coelitcrvente. 
tees wrand £ are spatial coefficients whieh determine the 
pPhase-delay and amplitude change. “they are products of 
rational functions of x and exponential terms involving 
Vekocd test mandswatio.oft depth of seurce to thickness. The 


expliciterelarions on their references for the three types 


of sources. gare as follows: 
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where Ry Rio S 3 and Sy are derived in the Appendix A. 

ihe displacements in tthe time “domain are obtained by 
an inverse Laplace transform following the methods of Cagniard 
(19390) Saekeris (1940), and’ de Hoop, (1960)... The detaits os 
the transtormation follow the method of Loneman (19621) “as 


Peneralivzed by Abramovici. GL97s). 
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Us ays The Weathered Layer Model 

The seismic weathered layer or zone is a surficial 
deposit: that has a profound influence on all studies. of the 
Ci1astic propagation of waves. Tie -propertres of the 
weathered layer are: 
(a) a low compressional wave velocity (300 to 2000 m/s) 
Ovid onOmarousity Low Shear velocktyeGloo0 to 1000 m/s): 
(pb). lie *enerey absorption ‘or Pow 0 CSOOUe: 
Ce)" "a Srererogeneous: composition wi th t hettrniteknes s! bviaranin ¢ 
Poms tO Lod meters. However the Lateral variation 
Weuciry tas 2a Long wave "Léneth “except at weoloericalhly 
Steurireart interfaces * 
Cry ea saree “Oharce in “acoustic impedence *atroses the inren— 
face separating the weathered layer from underlying 
consolidated medium. 
The propertires "arise “because "tire Troecks mMear "Ghe "surface are 
relatively urconrsolidated, mechanically “and. "chéemica ] by 
aatered and highly “aerated : In some ‘areas*'the base of the 
weathered tavyer ws "related tothe taverace™ level or the wares 
faple but even’ in’areas “where “the watten "fable iseate the 
Surrace there hs av’well¥defined*selsmie weathered Layer 
presumably caused by repeated hydration-dehydration, 
preciptration of Salts; fracturing and cecompostei on Ot 
minerals. 

THe effeet ofl fhe weathered *Vlayer*on a*strfatce’ retord- 
ing’ or’ on’ seismic waves “impinging on’ dt from greater depths 


Te" OULU wor ali proportion to its *rélative -thiekness or 
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PEreoviiweime cen Trowel at. The weathered layer produces 
anomalous scattering of energy, magnification of amplitudes, 


eicevaree enerey dissipation, particularly in the desired 


Pop frequency part of the Spectrum. Cr ectevas: (ay waite 
euade which converts!a large amount of the ‘energy into 
Sutrace waves and. generally introduces) undestrable rever-— 
Detaiton patterns into the seismograms. 

Mieraccurate determination of the weathered) layer ic of 
Vatiewitemeki ne static ‘Corrections for net lectron seismic 
Puucies and it has potential in the derivation’ of deconvolnu= 
Beene Operators ; The subject 12s of interest by ttecelt in. the 
Ser Ot Che water table, buried channel deposits and tire 
Sere tton oi coal and other suriiczial deposits whieh may 
Poet pect LO strip \minins. Because of the heterogeneous 
nature Of the weathered Layer itS. propertves may only be 
deLrermined by seismic Studtes employane Short source ‘to 
receiver distances. Bbasticc, Wave rTreeconrdings In the near 
Suoetinoermedbate idistance Tantvel are dadiicult 20 smtenpeet 
because setemolocists have been trained in and ose, tari eid 
jadeweerDor most Of Cher studies. In general, more care 
must be (taken in the type of instrumentation used for near 
Surface studies because of the diftitculties in identifying 
Severale kinds of phases. We shall «see that three componeurt 
recordings are essential since the radial direction often 
contains the strongest compressional or P wave arrivals. 


In fact. P velocity arrivals may be found on synthetic 


seismograms in the transverse direction in the near field. 
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Pepa eusde CO. Eh SEUdy (ol Seismace recoraanees, tor tue 
Bonvcture Of Shallow; layers, a Series ob \@mact synthetic 
seismograms have been computed using generalized ray theory. 
ime model to be used\is idealized in that it consists of 
Seoource within or on an elastic layer over a haki=s pace: 
Prewimcorpoeration of attenuation will peweonsigered in he 
mext chapter. Gegd approximations to va low velocity 
Peover over a half-space are found quite frequently in 
feurtmenecary basins. More complex layered Structure could 
besuodelled but tisually with objectionable approximations 
Pec al POrLthn or with prohibitively expensive computer 
operations. Before attempting these more elaborate models 
tie Considerable complexity in the results from these 
simpler examples should be explored. 

Poecheltollowine (synthetic Seismograms the complete 
Feri wom to oi ven ay toa selected fame arrem the fared 
re ray aly An important advantage of the generabized ray 
Method #£S, phat one may examine individual peneralized rays 
or Proups of them. Thiss-allows one to determine if the 
numerical integration procedures ane stable aud “accurate lox 


each generalized ray. Lt Ps «also pessip le Geo sdeternmare 


Whetweresesieniticant peak is due to a mingle ray ‘or, as 


Often heppens, because of a2 cluster of farrivatle. 
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Figure 1. Model of a layer over a half space with the 
source at a shallow depth. The model simulates 
wave propagation in a weathered layer. The 


source is a stress impulse with a triangular 


shape in time. 
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Masoure, 2 


Synthetic seismograms for a horizontally 
polarized torque as the source in the model of 
figure 1. The surface transverse displacement 
is shown in a three dimensional plot using an 
algorithm called ASPEX. The figure viewed for 
an azimuth of 340°>3 jan altitude of 91S" frome se 
horizontal at the center of the.graph to the 
observer and a distance of 25m from the graph 


to the opserver. 
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Sueveueat om Or at’ the surface. Tne ftirst ‘seismocramns 
(Fig. 2) are of the transverse GOmMpPoOneNnt: at distances of 
PO toro0 mm for a torque with the axis in the vertical 
airection . Bor -anvexperimentalist en SH torque is mot a 
werye interesting model since it is difficult or impossible 
BOMmDTOauce= in the: field. For the theovetieianeand procram= 
Mieteioeis the first sitep in producing a well tested 
Beagorithm since only horizontally polarized waves are 
produced. The results are of some lake ase Since apart 
PPoOnaverlocity scaling, an explosive source in a, liquid 
medium produces a similar seismogram on a pressure trans- 
Guceyr . Note the well developed head wave (labelled SS*= 
S,S58,) Separating strom the first rerkected wave. (ss). 
Phemoe ce. sounce Vorot far tereater “interes ec. tho the 
experimentarist. being avconcentrated vertical force 
produced by a jump ons normal stress at ta ‘particular depth 
Once: zo axis. A ‘vertical force ~us Gsed-=instedad, of va 
buried compressional pulse since a hammer impact, a 
MabpGeatony device or a. cylindrically shaped tube, of 
explosive chemical in’ a lightly tamped borehole produces 
Bidominantly vertical force. The difference in the displace— 
ment components for a surface source and one at a depth of 
5 m is shown in figure 3. The displacement for the surface 
source is similar to the classical Lamb (1904) solution up 
to the time of the reflected (PP) ray. Dramatic changes 


occur when the source is buried. The: amplitude or ithe 


Rayleigh wave is over one magnitude greater for a surface 
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Figure 3. Synthetic seismograms for a surface and a 
buried (5m) source for @ point vertical. torece 
The model parameters are given in figure l. 
The decay of the surface head wave sP* is shown 


to the right of the seismograms. 
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impulse. 

More importantly, there is a new diffracted pulse 
labelled sP* which travels to the surface as an S wave and 
along the surface as a P wave. This "surface head wave" 
was predicted theoretically by Petrashin (1959) who called 
it a shear head wave. Numerical seismograms for this phase 
which they called a diffracted S* were first shown by 
Abramovici and Gal-Ezer (1978). The Particle motion is 
Pinear and 180° out of phase with the direct Por S waves. 
Hiewsurtace head wave. does not occur if the source is purely 
compressional since no S waves are generated. From theore- 
tical considerations the surface head wave should decay as 
ae 2 SLOse VO Critical distance andvas Age at distances 
where the path in the lower interface is large. This is 
Contirmed by the modelling studies (Figure 3). 

For the buried source the emergent Rayleigh wave 
coecariyeprecedés and overlaps the direct shear wave... © bpis 
Dehavaoun can be anterpreted physically, as; due tov the 
Srrivableor the direct spherical Pswave front: at «the suntace 
near the source and generating a displacement which forms 
the early part of the Rayleigh wave. The interference 
pattern formed by the strong direct shear wave and the 
Raylei¢h waves will certainly create difficulties in inter- 
preting field recordings if only a single component is 
available. The interpretation of the displacement of 


generalized rays is facilitated 121 they are e2rouped into 


families defined by a parameter F 
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Figure 4, 


Classification (of 30sphystealirays inca 
generalized rays, having a unique value of u 
and Vv into the first three -bamities... vibe 
partial seismograms, using the vertical 
component at a distance’ of. 40m, and “ay verticas 
force buried at 5m, are used to illustrate each 


family of rays. 
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F is equal to the number of encounters that the curve makes 
with all available interfaces. In figure 4 all the rays in 
the first 3 families are illustrated together with their 
displacement for a buried source and a surface receiver at 
a distance of twice the layer thickness (40 m). The complete 
seismogram is the sum of all the families of rays up to any 
particular time. This may be viewed in figure 5 which shows 
the suite of seismograms from 10 to 60 m. As will be shown 
in the next section of the chapter, each individual general- 
ized ray has its component impulsive ray arrivals superimposed 
upon an “exponentially” growing tail which is not evident 
on the complete seismogram. The tail in each generalized 
ray may be thought of as due to Rayleigh waves generated as 
Ener spnerit cal) <and,:S wave ,frontsiinteract with. the surtace,, 
When all the generalized rays within each family are combined 
the tails cancel completely after the arrival of the last 
Taye "Aarelict, of the tails ts-phystcally observable by the 
non-zero offset of the displacement between individual rays. 
The partial seismograms in figure 4 reproduce faithfully the 
ground displacement including the non-zero plateaus between 
the impulsive arrivals. 

Note that the odd numbered families originate as up- 
going waves at the source and even numbered familites initiate 


as downgoing waves. The amplitude of phases in family 3 are 
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Figure 3. 


Synthetic seismograms for a point vertieal 
source buried at 5m for the model shown in 
figure ls. The three dinensioneal plot hacwmenve 
same parameters as figure 2 except that the 


reduction velocity, fs) 1000 "m/s: 
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= 1.0Km/s ) 


reduced time (a, 


Gievecoitently less’ than those of family 2. This as’ to be 
Parecred cince ‘the extra encounter is as a perfect reflection 
at the surface. In#¢field recordings this) will seldom occur 
ena sthe rays from family 3 would diminish by an amount 
determined by the roughness of the surface and dominant 
Boavelength of, wave and surface irregularities. Rays from 
Hagher order familiés~ have amplitudes much diminished. This 
momevident i1 figure 5 which also Paciudes rays irom <fami lies 
4 and 5. These arrivals are the ones not labelled with a 
code. 

Meo synithetic.ror a buried) vertical force (figure 5) 
one should note the prominence of the first two compressional 
waves on the radical component. These are the direct P and 
erewcouriace head, wave, -sP*. At intermediate, distances, 10 
poms Oometrom the Source, the shear (S) <arrivalsis more 
protanent ‘on the vertical receiver. There are a great 
mubciplicity of reflected and related head waves in the 


seismograms. These €xuct seismograms include. 20 |paixrsmort 


generalized waves originating as P or S waves from the source. 


The total computation time for all figure 5 seismograms on 
Ene Andant 470 Vi6 was 25 minutes at a nominal cost of S500. 
Unfortunately the integration noise increases greatly and 
the CPU time mounts rapidly if longer seismogram intervals, 
a@rter the first arrivals, are required’. 

Wienethe buried source 1s 4 Lorce in an arbitrary 


horieocateal direction all three components of motion are 
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Figure 6. Synthetic seismograms for a point horizontal 
source buried at 5m for the model shown in 
figure l. The ray paths) tor a-tew prineipan 


rays are shown on, the rienet. 
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excited. In the examples shown in figure 6 a unit impulse 

Sf Equal "strength im ithesx (=x) and y (=x,) dix ection. .was 
used. ~*ihe ‘source ‘acts likesasshearr dislocation with a stmess 
dascontinulty along °an open*surface.t Although) the sysotnem of 


equations is expanded conveniently into an SH set with 


transverse motion and a,P-SV+set with radial andi vertical 
motiomacboth setis fravelarrivals with Peand S»wave. velocities. 
Thus the transverse (SH) component shows an initial wave 
Prnivein? “wihtinraehtwavetvelocity atyneartand, intermediate 
distances (10-20 m). We have observed this with hammer 
seismograph sources in field experiments during studies of 
the weathered layer and the theoretical confirmation is 
ST a2eiryaine): 

Phe tdirect Peiwave isttheestrongest) marrivadewons ithe 
racial SCcomponent "while®théewshear® (S)earrayvall aise ther strongest 
On “bhe "vertical "component sr Excitationvoret hits tty ero 
Source is particularly) desirable’ since the RayVeiegiwaves 
are minor "and “the séismogram is dominated® by falrcilear tset of 
rerrected -tmpulses®S® A®Gcmparison* or figure: Si with awerti- 


Gal Vforce "and*fteure*6 with’a horizontal force *makescifene 


point *crearty. The surface head wave, sP*, is the only head 


wave with a well defined signature and it appears strongly 
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Zeo the cCrustal:-layer model 

Synthetic seismograms are computed for a horizontal 
point force in the 2 plane for the model shown in Fig. 
7. The model velocities and densities are typical for 4 
40 km thick crustal layer over a mantle half-space. 
Synthetics are computed for two different depths (figures 
cmeand 9) rand f0wodistances: of 20 to 60 km. “Inoalil. cases 
the receivers were placed on the surface on the x, ~axis. 
The x ,~component oflithe force ‘gives tise to radiata d,eand 
Venuicai, Ww, components of motion. The\-components of «the 
force along the X-axis gives only an angular displacement 
called the transverse component, V. As in refraction 
seismology the seismograms are plotted against reduced 
time (t-r/a,) where. ¢ fis othe actialUstimestne seconds. ote 
the distance from the source to receiver alone the x-axis 


L 


4 isthe velocity of compressional waves in the layer. 

Some of the rays are marked by the convention employed 
in earthquake seismology. An initial wave propagating 
toward the surface has the ray segment designated by a 
lowermicasesletter (por Ss). Avstiar is used “as (a-superccrtpe 
to indicate a head wave of that type. The solution is 
complete and includes all the generalized rays arriving 
within 23 seconds after the first compressional wave front. 
At most distances this included about 15 generalized rays 
starting at a compressional wave velocity and 12° startine lat 


a shear wave velocity. The integer yw is again a generalized 
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Hipurnes 7/5 


The crustal model of an elastic layer over a 


haliespace: The source is at 
1/2. Chi eure. S)§ orl Seite 


thickness. Velocities (a and 


while density (p) is in gm/mg. 


a depth of either 
9) of the dayer 


8) are in km/s 
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Figure 8. 


Synthetic seismograms for a source buried in 
the middle of the layer ot  ficure 72 tiem a. 
paths for a few of the principal generalized 
rays are shown on the right... The ficures ieee 
three dimensional plot for an azimuth of —340e 
an altitude of 18° from the horizontal at the 
center of the graph to the observer and a 
distance of 25 .km=eitvom the -oraph to ine 


observer. 
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ray index which may be interpreted as one less than the 
number of interactions with the available boundaries as an 
indicent P wave. Similarly v is a generalized ray index 
which specifies the number of boundary interactions minus 
one as an incident S wave. As in the previous section, 


one can define the integer F 


| OPER =r Ue Res oY ees aed . Poem wiles sieves lee 2 LD 


Piethe illustrations for fieures 8 and 9 the generalized 
ray parameters which needed to be included for the desired 
Svanmouulime owere 1h) s+ 60.0 cS 25,7. F os) 7. 

Although the source is due only to a horizontal force 
there is a prominent p or compressional wave generated by 
the first generalized ray (yv=0, v=0, F=1) at close distances. 
The vertical and transverse components become negligible 
wWithina distance equal to\ the thickness of the first layer 
but the radial component is substantial at all distances. 
This result is of significance when interpreting seismograms 
from horizontal sources in exploration seismology. The 
direct s or shear wave shows some interesting variations 
in amplitude on the radial component when the source is 
shallow (figure 9). It is also possible to see the head 


wave generated by shear waves incident on the boundaries, 


Ere one the cadial component. 
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Figure 9. Synthetic seismograms similar to those shown in 
figure 8 but withthe horizontal pointe source 
at a depth of 1/5 of the Layer thickness etme 


three dimensional plot uses the same parameters 


as in figure 8. 
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Figure 10. Phase relationship: in the near “field tore 
buried horizontal source in a, layer) overwed 
half space. The transverse (V) radial GQ) @and 
vertical (W) component are shown for a receiver 
at a distance of 20 km. The source is buried 
at a depth of 8 km. The model is showng an 


fi gure. 
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For horizontal force sources the vertical and radial 
components are close to 0° or 180° in phase for most 
generalized rays at medium recording distances. However 
surface recordings close to the source show large phase 
biflerences at times close to the direct s wave arrival. 
This is shown more clearly in figure 10 which is an 
enlargement of the seismogram at 20 km in figure 9. The 
radial and vertical components following the s phase show 
that a small "Rayleigh" wave has been generated from 
COMtimeDuUtLTons Of inteprals'’ with indices )= O;6v = 0sCr = 1). 

Since the final seismograms were computed from a 
linear superposition of generalized rays, it is possible 
to disassemble the contributions. Examination of the 
individual generalized rays allows one to make a better 
physical interpretation and to check on the numerical 
stability when evaluating complex integrals on a digital 
computer. In addition it may be possible to compute an ad 
hoc synthetic seismogram including attenuation by treating 
each generalized ray individually and incorporating atten- 
uation in a process involving a discrete Fourier transform 
Cohapeer 3). 

The generalized rays may be grouped using Eq. C221) 
according to families based on the value of F. Thus, F = 1 
defines the first family which include the direct p and s 
arrivals and their associated Rayleigh wave due to the 


interaction with the free surface. Figure 11 shows the 
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Figure 11. Synthetic seismograms decomposed into families 
of generalized rays for a distance of 40 km 
and a point horizontal ‘force in the middtemo- 
the crustal layer. ‘P and 5S’ indicate thaterene 
ray starts as a compressional or shear wave 
from the source. All the partial seismograms 
start at a real time of -6.8s.- Note that whe 
all the generalized ray contributions in eager 
family are summed (2), the exponentially 


growing tails cancel. 
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individual generalized rays for each component of motion 
and their summation into a family of rays. Each individual 
generalized ray shows the principal arrival as a P or an 

S wave but it also shows a tail which does not vanish with 
Pane. (ihe tails are related to: the generation of an inter 
face (Rayleigh) wave as the curved wave front expands to 
eoversa COntinuoussrange of’ ray parameters. As shown..in 
figure l1lthe tails disappear when all the generalized rays 
fromecach family (F =F v + 1.= constant) have arrived 

and are incorporated into the seismogram. PON saceen@ that 
from a physical point of view F is a measure of the number 
of encounters each ray has had with the available boundaries. 
The classification of generalized rays into families 
appears to be useful and should allow one to form an alter- 
mate method of truncating a seismogram when it is too 
expensive to incorporate computationally all the rays to a 
particular time. 

The synthetic seismograms should be of use in inter- 
preting field results close to small earthquakes where the 
source can be modelled by a horizontal stress discontinuity 
at some depth within the crust. The importance of the 
direct compressional wave and many converted phases has been 
illustrated. The results should be of interest in explora- 
tion studies where a shear source is generated by a hammer 
mechanism or by a vibrator. Since the solution is analytic 


and the numerical results are complete up to any desired 
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time the results are useful in checking other methods 
which may be adapted more easily to a complex elastic 


layered medium. 
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CHAPTER 3 


leet roduction 

It is becoming increasingly important to include the 
effects of anelasticity in the computation of synthetic 
seismograms and in other seismic modelling techniques. A 
method of introducing the effects of attenuation and 
dispersion into the synthetics will be shown. The general- 
ized ray theory is used to calculate exact synthetic 
seismograms for a layered homogeneous isotropic medium and 
tiewermects Of anelasticity such as attenuation and dispersion 
of seismic waves are modelled using the linear theory of 
viscoelasticity combined with Futterman's model (1962). 

According to the generalized ray theory the response 
of a layered medium to a disturbance is expressed as a 
superposition of individual generalized rays. By decomposing 
the seismogram into rays, one can study each ray individually 
and the effects of absorption and. dispersion can be incor-= 
porated in the frequency domain using Futterman's theory. 
Eneaddition, the effects of viscoelastic intertaces ane 
Faken) into account by calculating reflection coefficients 
for anelastic media. Therefore, by summing all the rays a 
synthetic seismogram which includes the effects of anelas- 
ticity can be computed. 

Futterman's theory provides an excellent model of 
attenuation and dispersion that is in good agreement with 


experimental data (Wuenschel, 1965). It is based on the 
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principle of superposition and on the linearity of 
absorption coefficient as a function of frequency. Reflec- 
tion and transmission coefficients for anelastic media are 
calculated using the correspondence principle. Lockett 
(1962) described the method for setting up the equations 
which must be solved for the reflection and transmission 
coefficients in the anelastic problem. The boundary 
conditions are the same as for the elastic case and the 
only change is that elastic modulus, phase velocity and 
wave number are complex functions of frequency. Buchen 
Gihoa)ttand! Borcherdt (19)73)° set Aip “the *ceneral™ theoretical 
framework for plane waves in linear viscoelastic media by 
considering inhomogeneous waves, where the direction of 
propagat lon Wis difterent* fron ‘the ‘direction of “attenuation. 
Schoenberg (1971) solves the problem of a plane harmonic 
wave impinging on a plane interface between an elastic and 
a linearly viscoelastic medium, Borcherdt (1977) gives a 
complete theoretical description of SH-waves in anelastic 
media and derives reflection and transmission coefficients 
for the SH-waves. Krebes and Hron (1980a,b), following 
Borcherdt's formulation, calculate reflection and trans- 
mission coefficients for the SH-anelastic case and compare 
them with the coefficients for the perfectly elastic case. 
They use these coefficients for computing synthetic seismo- 
grams for SH-waves in a layered anelastic medium using the 


asymptotic ray theory. In this chapter we consider 


attenuation for both the SH and P-SV cases by a modification 
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of an exact and complete solution obtained from generalized 


ray theory. 


3.2 Attenuation and dispersion along the path 

A complete description of Futterman's model of atten- 
PattOonwasvei Ven. imsthe wZirst chapten. skauations (1.31), 
CL. 2) wand. (1.33) describe the attenuation model used in 
Citic thesis. (They provide the absorption coefficient, the 
phase velocity and the dimensionless factor Q as functions 
Cimipmequency. “Then, the effects of attenuation and 
dispersion along the ray path can be easily incorporated 
into the frequency domain by multiplying the amplitude 
spectrum of each ray by an exponential decay and by adding 
2-pnase shirt in “the phase “spectrum? "The “nay “ls ‘obtained 
in the time domain using an inverse Fourier transform 


CEq #*2.29)*5 


S23 4tihe effect lof “the viscoelastic) interiace 


The stress-strain relation for a homogeneous isotropic 


linear viscoelastic medium is given by 
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where ows istthe®stress tensor; 7A and uw are the complex 
ij 
Lamé parameters and e,, is the stress tensor (Borcherdt, 
1j 


EOF 3 P1977 )4 Substitution of (3.1) into the equation 


of motion gives 
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[aCt)+u(e) ]*a(viva) + w(t)*d(y7a) = pf (3.2) 


after some calculation where the symbol * denotes the 
Stieltjes convolution. Taking the Fourier transform of 
(3.2) and writing the transform of the displacement vector 
in terms of the Helmholtz potentials; we obtain the well 


known Helmholtz equation 


VF + KF = O C33) 


where K is now complex. 

Equation (3.3) is similar in form to ‘the corresponding 
@quation of linear elasticity theory, except that the 
Lamé constants are now complex. Therefore, any formal 
solution of Helmholtz equation in the theory of elasticity 
offers a corresponding solution for a linear viscoelastic 
DOdy~eiltechne elastic) moduli that 7occur jinethe elastic 
solution are replaced by the corresponding complex modulii 
(Correspondence-Principle; Ben-Menahem, Singh, 1981). Thus, 
reflection/transmission/coefficients for a viscoelastic 
interface may be obtained from the analogue in Che elastic 
case. In appendix B we show how complex modulii are 
calculated using the phase velocity and the attenuation 
factor as functions of frequency. 

Figure 12 shows the amplitude and phase for the PP 
and PS reflection coefficients. The solid lines represent 


the elastic case while the dashed lines correspond to the 
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The amplitude and phase for PP and PS plane 

wave reflection coefficients are plotted against 
the angle of incidence 6. The solid lines 
correspond to the elastic case while the 

dashed lines correspond to the anelastic case. 
The velocities, densities and Q values are 
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viscoelastic case. In order to examine the frequency 
dependence of the reflection and transmission coefficients 
for a viscoelastic interface, we have studied several simple 
models. We found that these coefficients are not sensitive 
functions of frequency unless the values of Q are very low. 
Therefore, the effect of a viscoelastic interface can be 
modelled by calculating frequency independent reflection 
and transmission coefficients for anelastic media. As 
Peeuneni2 sindicates, in order to take into account. a visco— 
elastic interface, we have to apply a correction to the 
amplitude spectrum of the pulse and to add another term in 


the phase spectrum. 


Bo lhe attenuation algorithm 

Procure 15-16 1a flow) .chart wiien vilustratesatie 
algorithm used in order to introduce the effects of atten- 
uation and dispersion in the seismograms. The synthetic 
is first decomposed into individual generalized rays, and 
each ray is transformed into the frequency domain using a 
fast Fourier transform (FFT). The absorption and disper- 
sionsatone the ray path is then introduced using 
Futterman's model. The effect of the viscoelastic interface 
is also taken into account by calculating reflection 
coefficients for anelastic media. 

The procedure for one generalized ray is illustrated 
in Figure 14 for the case of a horizontally polarized head 


wave and reflection near the critical distance. A Fourier 
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Figure 13. Flow chart showing the algorithm used in order 
to introduce the effects of anelasticity nee 


each generalized ray. 
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Figure 14. 


Steps involved for the inclusion of attenuation 
into one generalized ray. The horizontal 
distance between source and receiver is 30m. 
The source is an SH-torque for the model of 
figure 15. (a) Partial séismogram tor ene 
generalized ray (v=1) in the case of a layered 
elastic medium. (b) The amplitude spectrum for 
(a) as obtained from ia) -fast Fourier transtorm 
algorithm.  (Ce)) The phase? spectrumyfor aca 
(d) ‘The unwrapped hase for Ce). <e) iP ameueen 
seismogram with the effects of anelasticity. 
(£9) The amplitude spectrum for the attenuated 
pulse. (g) The velocity dispersion curve for 
the first anelastic waver an tieure. 15). 

(h) The unwrapped phase for the attenuated 


pulse. 


84 


= 


= 


(4) 


SO- 
it 
CC 
(e) un 
SO 


AN 


——— 


ae 


LL 


~ 35 
24S) 


. 
\ 


15g 


\ 


33 


‘tHure ip: 
gunerolided 
ei@atic ned 
: ha 3 | 

iad ae obtatoae’d ; 
ei cweatios Soe). he 
- i’ = s 

Len i hese fier (4)~ fe) Baw 
y@eAgi a aiek the eee ts ot @nelaut eae 
ite a a my Dee 


Pies | 
a 
= a 
~ ‘ 


PEs 
_ 


SOURCE. - 


” 
@ 


a? 7 
Ppl: iol 


if 
ae 


% 


MODEL 


TL ee 
a 2 | is 
re : as 
i vf : ah r} 
a on ee t 
i 


viebodt!) tay 
SD aeiugavts) a dddw « 


Po Hetaotbak 2 mots 


a 
ie 
fi 
y) 
l 1] 
j 
! ww 
ste 
= 
a 
- won ie) 
_ 
oe 
os 
Tomy 
tn, 
i “=a 
ronan 
u a ae tm 


» o > rr oe eal Sie nr 


i 


hy 


‘ 5 rey 
i @ ' Sry 


re in 
ah Sh3N 08 


pe. 
i Ma ed 
Maem! 


re cust 
: a 
C4 
iB! 
i 
‘ 
ia | 
fs 


2%) . ~ 


i" cf) eave it 
~~ 


Figure 15. The anelastic weathered layer model. The 
source is a stress impulse with a triangular 
shape in time. The attenuation is indicated 


by thes quality efactoreu. 
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transform is taken of the original seismogram for one 
generalized ray. Its phase is unwrapped and then a phase 
term is added using the phase velocity c(w) in equation 
(1.32). The modulus or amplitude function is multiplied 
by an exponential decay involving the absorption coeffi- 
cient. The effective change due to the viscoelastic 
feerect+on coeitacient is introduced similarly: and ‘the 


final pulse is obtained by an inverse fast Fourier transform. 


3.5 The weathered layer model 

A model which simulates the weathered layer is shown 
foo tagure "15. It consists of a single homogeneous isotropic 
layer overlying a homogeneous isotropic half space. The 
compressional (a) and shear (8) velocities along with the 
densi tres Go) sanditche 0 values ‘are aiso indicated. 

Synthetic seismograms for the elastic and anelastic case 
are, calculated for two different sources, an SH-torque: and 
a vertical stress discontinuity. 

Figure 16 shows the synthetic seismograms with and 
without anelasticity for the weathered layer model using a 
point SH torque as a source. The displacements are presented 
in a three-dimensional graph with one of the dimensions 
being reduced time. The seismograms are at intervals of 
5m for a surface receiver at distances of 10 to 60 m. The 
direct wave S, the reflection SS with the corresponding 


head wave SS* (=8,S,8,); and the multiple, sSS, are shown. 
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Figure 16. Synthetic seismograms for an SH-torque source 
for the elastic and anelastic case for the 
model of figure 15: “Site aleorithm AsPExX ic 
used to produce a three dimensional plot for an 
azimuth, of 340°. tan Sei tude of (hs - et ronmene 
horizontal at the center of the graph to the 
observer and a distance of 25m from the graph 


to sthe observer. 
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Notice the decrease in amplitude and the widening of the 
pulses in the anelastic version. The decrease in ampli- 
tude is due to the exponential decay while the increased 
pulse breadth is due to the loss of high frequencies and 
the svefliects,of dispersion. 

ihewe-oviecace “is TViustrated: in ‘hieures,) (7. 1 smand 
iJon which aa burbed point vertical force Js uged ac a 
source. Figure 17 shows the vertical and radial components 
of individual rays resulting from the decomposition of a 
seismogram from a vertical stress discontinuity. The 
horizontal distance between source and receiver is 20m. 
The rays of the first two-families (F=1, F=2) are shown. 
Each family is completely defined by an integer F which 
indiacates -the number of interactions. of the ray- withthe 
two boundaries. A Fourier transform of a tapered: ‘form ‘of 
each generalized ray is obtained and the attenuated pulse 
te-obtained following the procedures outlined in Figures 
13%and’ 145° The complete synthetics for the vertical and 
radial component are shown in Figures 18 and 19 respectively. 
The seismograms for the elastic and anelastic cases are 
plotted from 10 to 30m at an interval of 4j73m. Por; short 
distances between the source and receiver the path of the 
head waves was taken to be the same as for the reflected 
waves. Similarly, the generalized ray which includes the 
direct wave and the Rayleigh wave was treated identically 


for the purpose of attenuation. Note the dominance of the 
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The vertical and radial components for the 
generalized rays making up the first two 

famidies CF=1y7 h=2)e CP)esande (S)pandt cate 
beginning from the source as compressional 
shear waves respectively. The sums of the 
generalized rays making up each family are 
labelled, by a-°2..1 Noteuthe ‘cancel lationset 
tails in the partial seismograms for each 


family. 


rays 


and 


the 


= Ss <20ms— 


92 


VERTICAL 


RADIAL 


$=4 7 
= Aa, 
. 
. Sc 
oo ee ee Nee ee ffi) — we 
Ce > fsa “ 


kee} ee ee +P) “ah — horace 
O=« WJ [= a -t 4 ; 
env sical ase tndted compan ay ' 


vere 2 | 
- ya oars <comaiiacd ties eantee ap She cirae! 


i= y Oy 


Te ay ohition (¥=t, Ped). 68) ee te) 


i 


te cae errrer . tebe “the saurce ae ae KER: 
ee ay cal = : : 

ae oy wt weveg eee Sve ty. The @ 
“i £@ ie PAKIng up need fat ly @ 


~ 


Pl (ied tm 4k. Mbte tie capcelletion ; 


paretal So \seegeaems for 
lst “wermmnceme 
tw Oe ae 


. +e 
ay - 


a 
a 8 


— 7 7 
ae We 
at ft ? 
y | 
: : 1 ' 
e ;f j 
> 4 
i 
f 
‘ ‘ 
i = aneiasric 
+" o 
— Sa ee - eae * 
o q es — a a. ¢ 
a = rae. _ | tt oN 
7 S ie, a 
a tees ~ Mie RP Atieers » nt 
¢ ied . SRP ear as 
” _ ee ae . \ i eT 
—~, a *, we ee a . 
Th hae ll ae ie ae — *, ~ : a — a. a 
0o—..... ~— ee ‘ 7 at io S = im - — 
wa sip - ee ee 
i od any on { —_ i a ehenial —— — 
pe ar eae | \ ite “ _ —o ~ ate - a. 
Ad cl , = ~~ ? 7 oe oe ah eng 
aA =, >» ane _ — 4a ee 
———— "| a aes - ‘- io el ieee q 
; = —, , > oe *- tt cat anes, 4 
meee |) ee | 
: ¢ J 4 ; 7 
' - 7 oor 
i t a "1 
; J i i » ; 
b ' : 
/ | ' 
. efdOGTE 302 caancaron (astotev sat 28E aaygtt 
gnitog tacisrey's 103 raneaqras stitev edt <6L 3 
oF) 
Se ? *; aid * 
wee Ssee cided (sta -tue-risculs.. 5.70) p22 1 Oaeg ws $$ 
= j a . y _ rr. 
. ~ =e J on if Nantes 
3 a the - /‘ abe ta awk 2 —, 
Wee. tl wriigti 30 latte iay! PI S¢ 3496 qt me Nii 
; = —_ ee ee Sey eed \ 
4 e- ON ee 
idistonce et ee . a © os ee ee a a | a 
; : j — 2 a ; , sees oo 
Of-s200)-4>c4..479 5 Sede Ta PRs. O40 /E259 
e f Ping sing ree “a Pd we ee _ , ier Oe 
y P — 
NN pre ay a rr le ne , oa 
wy ‘ammenities bi 7 ~ ee ed ee ee i al ere 
ie ame O10 \ ; P a ae, ae 
i = —— WW 
A ‘ \ — oe | nl il 7s ee A A 
. j= ae “ i" wil — - ~~ 
Vv 7 = mm -_ ae sail 
‘ i a _ a a 
; ra , — — ere... -— = 
’ al — - “a 
ra ae 73 —ert ad a 
‘ee é0n 
| 
{ 
ods Heren, Of oe LOR | 
hs 
Zo ap 
} 


Figure 18. The vertical component ior jay vertical apoiae 
force for the elastic and anelastic case for 
the weathered layer model of figure 15. ASPEX 


uses the same parameters as figure 16. 
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Figure 19. The radial component for a verticals polnteaenm. 
for the elastic and anelastic case for the 
weathered layer model of figure 15. ASPEX uses 


the same parameters as figure 16. 
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Rayleigh wave as the high frequency components in the head 
wave are strongly attenuated. The head waves and 
reflections are strongly attenuated and show pulse widen- 
ing due to the filtering of the high frequencies and the 


effects of dispersion. 
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CHAPTER 4 


4.1 The "Alberta" Model 

Many sedimentary basins may be modelled, in a first 
approximation, as a thick, relatively uniform, low velocity 
layer overlying a thick high velocity section. A good 
example is the Alberta basin which extends into the Rocky 
Mountains on the west. In the foothills it consists of a 
3 km section of alternating shales aad sande of Mesozoic 
age with a mean velocity of 4200 m/s (13800 ft/s) overlying 
Paleozoic and Precambrian beds with a mean velocity of 
ee 00em7 sce (21000%f c/s) .9 The high velocity in the Paleozoic 
rocks is due to a large proportion of limestone and dolomite 
while the Precambrian section, of similar velocity, consists 
of 2 billion year old gneisses, metasediments and igneous 
avtLMostVess Chortert et ales 1962)" 

The mean velocity contrast between the two sections 
is) so0 large that special techniques may be used in mapping 
the contrast as discussed by Blundun (1959) and Richards 
(1960). Nevertheless, the identification of the phases used 
in mapping has been difficult because of the interference 
of many types of wide angle reflections and head waves. The 
computation of exact synthetic seismograms is also very 
difficult and has only been accomplished recently with large 
digital computers. Our own approach to this problem 
involves generalized ray theory and a Cagniard-Pekeris 


inversion because the method yields the complete and exact 
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vertical point force as source situated at 
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solution, and also allows one to decompose the results into 
individual generalized rays for detailed analysis. 

A concentrated vertical force at a particular depth 
is used instead of a buried compressional pulse since a 
vibratory device or a cylindrically shaped tube of explo- 
sive chemical in a lightly tamped borehole produces a 
dominantly vertical force. Mathematically it is simulated 
by a jump in normal stress across a boundary at the source 
bevel placed on the z axis. The model of a buried source 
in a layer over a half space is a gross simplification of 
the actual crustal section in the Alberta basin. A more 
complex layered structure could be modelled but with 
objectionable approximations in the algorithms or with 
prohibitively expensive computer operation. Before 
attempting these more elaborate models the complexity of 


the results from these simpler examples will be explored. 


4.2 Synthetic seismograms 
The basic model used is given in Figure 20. Other 
velocities will be illustrated after the results for this 


case are shown. To obtain a perspective on, the relative 


amplitudes exact seismograms were computed for distances 

ef & to 8 km from the source (Figure 21). The Rayleigh 
wavee dominate the recording and the earlier portion of the 
seismogram with the head and reflected waves is not seen 
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Figure 21. Synthetic seismograms for the model shown in 
figure 20. The vertical and radial componvene. 
are superimposed to show the phase relation for 


distances of 4 to 8 km from sthe sounce. 
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Synthetic seismograms for the model of figure 
20. The head wave’ from the half space” isea 
first arrival beyond a distance of’ 13) km. une 
algorithm ASPEX is used for the three dimen- 
sional plot withwan azimuth, Of 10340) .eean 
altitude of 18° from thewhorizontal 4 eethe 
center of the graph to the observer and a 
distance of -25Fkm frome the -25eph to jtue 


observer. 
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but at critical distances, we will restrict the display to 
the time before the surface waves arrive. Figure? 22.shows 
the vertical and radial displacement for a vertical point 
HoLrcegat distances’ of 6 to: 20 km (5 to 12.4 miles). These 
may be compared to the field recordings of Richards (1960) 
in the foothills belt of Stolberg over a distance range of 
Bo eO oi. O° kin, 

The major event on the seismograms in Figure 22 is 
Wiemrmeretlection, from: the half space. ~1t <tseioilowed 
immediately by two generalized pPP and sPP, which also have 
strong amplitudes. These two rays reflect off the surface 
Glose tio thes source. Their ray. paths may be seen in 
PLeure 4 . Their visible effect depends upon the depth of 
the source and the roughness of the reflecting surface 
which could be modelled as a parameter. The theoretical 
relative strength of the interfering, phases may be seen 
an PLeure 2:3.. 

The synthetic seismograms in Figure 22 also show 
strong arrivals for the SP and PS phases (see the second 
family, F=2 in Figure 4). The SP ray occurs only at close 
distances because a critical angle is reached when the 
compressional wave is travelling horizontally. At distances 
of 18 km the PS reflected ray may be interpreted. Note 
that the pSS* head wave in Figure 22 may be well enough 
detined to be of value in structuralemodelling. These head 
waves belong to the second and third families with ray paths 


labelled pS,P,8y and Si Po5y in Figure 4. 
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Figure 23. The main reflected and head waves for the 
model in figure 20. Three generalized rays 
are shown (v=1, v=0 (CP) pw=25, v=0 CP) 2 eze 
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Figure 24. Synthetic seismograms for a vertical pointe. onrce 
and variations on the model shown in figure 
20. The only parameters that are changed are 
the P.and’S velocities “in the halt space <a se 
P velocities vary from 5600 te 6800 m/s and 
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The effect of changes in the velocity of the Paleo— 
zoic half space is shown in Figure 24. The interference 
of the direct P wave, the PiPoPy head wave, and the PP 
reflection makes it imperative that an exploration program 
be carefully designed based on the thickness and velocities 
encountered. Note that the radial component has well 
defined impulsive phases and it would be advantageous to 
record on horizontal motion seismometers when carrying out 
seismic surveys for wide angle reflections. The results 
from Figures 23 and 24 show that one should keep the source 
at as uniform and shallow a depth as possible to avoid the 
gnterterence created by sPP and pPP rays with the primary 


wetlection, PP. A pattern of multiple sources would also 


help in reducing the interference of converted and direct 


shear waves. 
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CHAPTER 5 


Conclusions 

It is possible to obtain exact Synthetic seismograms 
for various kinds of impulsive forces using generalized 
ray theory and an algorithm that performs an inverse 
Laplace transform. Decomposition into generalized rays 
allows one to analyse component phases and interpret the 
results physically. In addition, the numerical Stapa ity 
of the integration for each generalized ray may be monitored 
independently. Synthesis of the generalized rays into 
Families is also of value in. the study of groups of srays 
and the influence of each additional interaction with an 
interface. 

The individual generalized rays are examined and 
modified and the effects of attenuation and dispersion can 
be incorporated in the frequency domain. Futterman's 
theory is used to modél the anelasticity along the ray path. 
Tae effect of the viscoelastic jinterface 2s also (taken into 
account by calculating reflection coefficients for anelastic 
media. Novel to this thesis is the incorporation of 
attenuation along the path and at the viscoelastic interface 
for the .P-SV case. In addition, use of a Laplace transform 
in the formulation of the problem insures that the final 
synthetic seismograms have a causal source pulse. 

In the near-field and at intermediate distances the 


field seismograms may be difficult to interpret without the 
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use of model results. 


trated with a buried source in a 


AKAs} 


The weathered layer case is illus- 


layer over a4 half space. 


A new surface head wave (sP*) which is not generated when 


the source is compressive is a prominent second arrivals, 


particularly on the radial component. 


a compressional wave velocity is 


A wave travelling with 


also Visible on the, trans— 


verse (SH) component for an impulsive horizontal force in 


areerDitrary direction. At near 
the direct compressional wave is 
component while the direct shear 
vertical seismometer. 

The same algorithm may be 
lower velocity layer over a high 


model), 


a solution seen in some sedimentary basins. 


and intermediate distances 
much stronger on the radial 


wave is best seen on the 


used to study a thick 
velocity section (Alberta 


The 


complexity of head waves and reflected arrivals near the 


critical distance makes it imperative that model studies 


accompany the interpretation of data recorded at wide angles 


of incidence. 


Most of the synthetic seismograms are plotted in a 


three-dimensional way using an algorithm called ASPEX. 


Seismograms for different distances are plotted against 


reduced time and the main arrivals can be easily identified. 


We found that this particular presentation of seismic data 


is very helpful and it may be of use in seismic studies. 


We also found that the algorithm ASPEX gives a better 


graphical representation than the well known algorithm 
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The models that we have studied are rather simple, 
in that they consist of a single layer over a half space. 
The obtained exact seismograms show the complexity of the 
results. More complex layered structures could be modelled 
but usually with objectionable approximations in the 
algorithm or with prohibitively expensive computer opera- 
tions. Before attempting these more elaborate models this 
thesis examined and studied in detail the considerable 


complexity in the results from simple models. 
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APPENDIX A 


Mie ene rTOdUuction 

The displacement components for a concentrated 
horizontal force in an elastic homogeneous layer on top of 
an elastic half space are given analytically in terms of 
generalized rays. Before going into the mathematical 
details let us give an overall picture of how the solution 
is obtained by reviewing briefly the various stages 
involved. 

The starting point is the representation of the 
solution for a stress discontinuity along a given finite 
surface © in a homogenous elastic solid in a convenient way. 
The solution in the frequency domain was given by Maruyama 
CI963) "for a discontinuity in: bothestress and displacement, 
in terms of solutions of the wave equation that are spheri- 
early symmetric with respect to points=on 2. Thus, the 
displacement components are written as finite combinations 
of integrals of such functions over =. When the medium has 
horizontal boundaries, the spherical symmetry is replaced 
by cylindrical symmetry with respect to the vertical through 
points of = and the boundary conditions are satisfied 
accordingly after rewriting the solution. This task was 
accomplished by Ben-Menahem and Singh (1968) and Singh, 
Ben-Menahem and Vered (1973) using the so-called Hansen 
solution (Stratton, 1941). We followed closely the treat- 


Ment used by Ben-Menahem and Singh in their 1969 paper and 
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split the corresponding non-homogeneous system of equations 
of motion into one involving P-SV waves and another one 
representing SH-waves only. These solutions were expanded 
apeseries Of generalized rays, each ray being inverted in 
the time domain by a generalization of the Cagniard-Pekeris 


method (Abramovici, 1978). 


Azw Lie formal ‘solution 

Consider an elastic solid consisting of a homogeneous 
layer of depth H over a homogeneous half-space (Fig. 1) and 
assume that inside the layer there is a time-dependent 


stress discontinuity of components 
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along a finite open surface 2. 


Our problem is to find the displacement vector 
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2. vanishing of stress on the upper surface, 


eg=aOutor z5=G0 (A.4) 


6— Bcontinuity of displacement and stress at 


Pipe rbace., or layers! 1 and 2 


U, = Us and t1 = Zo for z= H CA 25) 
a pthe! radiation condition 

useu0iferiz 2a (A.6) 
ie the source condition: when approaching £ the 


displacement ui tends tos, the displacement 
corresponding to the given stress discontinuity 
in a homogeneous medium with the same density 


and Lamé parameters as the layer. 
The Laplace transform of thetdisplacement 
u = u(p) = | u(t)e Prat CAT) 
is found as a superposition between the transform of the 


solution corresponding to the source s(p) and the transform 


U(p) of ‘the solution of (CA.2) chosen so that 
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u=-U+s (A.8) 


Saticties the transform jof conditions (A:4)—(A.6). The 
initial conditions are taken care of when finding U as a 


solution of 
Z 2.- - 
(uV"-pp’)U+(Atu) V(V-B) = 0 (A.9) 


and the source condition is met due to the superposition 
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A2a. The source solution in cartesian coordinates 


According to Maruyama (1963), the transform of the 
displacement components for a stress discontinuity along 2 


in a homogeneous medium are 
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zZ 2 
Here ho de hs ) are spherical Hankel functions. of the 
second kind, oe is the Kronecker symbol, and 
5 
Kee Dip es pos, one epee ere ou Re =e CRG Rua) CAR) 
J J a) Pe 


0,8 being the P and S velocities respectively, ze being 

the coordinates of the receiver and the summation being 
used whenever two indices in the same term coincide. 
Maruyama gave the Fourier transform of the displacement 
components so that in order to obtain the Laplace transform 


we change iw into p. 


A2b. Free solutions of the momentum equations. 

in ‘order to be abies to satisftfy the Doundary conditions, 
One must look for vector solutions: of) (A.9) that are 
separated in cylindrical coordinates. Guided by Hansen's 
solutions for the Maxwell equations (Stratton 1941), Ben 
Menahem and Singh (1968) considered the following indepen- 
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and a the Bessel function of order nm. 
These vector solutions represent upgoing and downgoing 
+ o 
eviindrical waves; La corresponding to P-waves Mo to SH- 


ic 
waves and Na to SV-waves. The general solution of (A.9) 


mreeerms OL Such waves is therefore 
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Ther tactors h.,kifand. © were introduced here for convenience, 
ai etc. being arbitrary functdonsrot i. 

It seems that all we have to do now is to write the 
source term in the same form, to add it to U and impose the 
boundary conditions. It may be, however, not soconvenient 
to proceed in this manner, even if the source term will 
turn out to contain only a few upgoing and downgoing waves 
and therefore only a finite number of terms will appear in 


the sum s0CA..L 5). The boundary conditions will result in 
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vector srelations between L ; etc. at z = 0 and Zu= H 


B I+ 


involving six coefficients in the layer and three in the 
half-space. These coefficients depend only upon ¢ and it 
may not be so easy to eliminate r and $6, the source of the 
trouble seemingly being the fact that although Li, ever, 
are independent solutions of (A.9), at each point they are 
not linearly independent, their number being more than 
three. 

Ben-Menahem and Singh showed an elegant way out by 
expressing the solution in terms of the following system 


of independent vectors: 
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The upgoing and downgoing waves are expressed in terms 


of this system as follows: 
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fo being the derivative of ee With irespect to 2. “The 
j j 
stress on a horizontal plane is expressed in a similar 


manner as a sum over m with: 
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A®e.. The source term, in cylindrical ycoordinates. 


Following closely Ben-Menahem and Singh, the 
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components a given by "(ASL ane 4first expressed in 


terms of the spherical wave functions 


mye = bay (-ipr/c) P) (cose) cosj¢ GAS23) 
9 : 
ae = eed (-ipr/c) p> (cosé) sinj¢ 


where 9,8 are spherical coordinates centered at the running 
point (E1268) on + and 2 are the associated Legendre 
functions. Using the Erdely identity expressing the product 


yor as an integral of scalar cylindrical wave functions, 


the components Sok are obtained as follows: 
7L = | a(t)J,(or)dt + cos2¢ | b(o)J,Cor)de 
0 0 
84 = sin2¢ | b(t)I, (or) de 
0 
S547 © cosd¢ | e(t)J5,(or)dt 
i (A.24) 
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d being the depth of the source. 

The expression for the displacement vector corres- 
ponding to the source in terms of the harmonic vectors 
Ba? Gs Pa is obtained using also the relations giving the 
connection between cartesian and cylindrical coordinates as 


well as ‘Some relations between Bessel functions. The 


Result is 
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Thus, the displacement due to the source is expressed 
as an integral over 1 of a sum of the form (A.15) including 
Docyetwor terms: | m =O and m = 1. If the surface ‘2 is 
reduced” to a “point, ite: the source is a’ space™ concentrated 


Boncesathe integration with respect to 2 will result in 


A SOE representing the strength of such a concentrated 


force in the directions of the coordinate axes and the 


Surcace, integralgin (A.26) will’ disappear: 
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The expression of s(t) in terms of L etc. was given 
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by Singh, Ben-Menachem and Vered (1973). 
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free solution represented by (A.19) it is convenient to 
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For a horizontal force in x,-direction with m 
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15 = jo = Ji = 0 
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For a horizontal force in the x.-direction 
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Thestotal.displacement-forea point forcessobtained 
superposing the source and the free solution is therefore of 
the following form (Ben-Menahem and Singh, 1968): 
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ce 


in the half-space: u,(P) = Re | u,(t)ede 
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A528, being respectively the P and S velocities in the 
half space. The coefficients arr 8 are obtained from 
the boundary conditions at. z = 0 and z = H,; taking into 
account that the vectors BE’ Pn? oe are linearly indepen- 


dent. The linear system obtained splits into a third-order 
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System for the coefficients cue 


and representin 
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the SH-wave, and a sixth-order system for the other 


coefficients, representing the P-SV-waves. 


A3. Ray expansion 
a. Pure Shear Waves 

Combining the source with the free solution, solving 
the linear system determined by the boundary conditions, 
expanding the solution in a geometric series, and grouping 
conveniently the terms of this expansion, we get the m-th 


term of the displacement at the surface in the form: 
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where By is the S-velocity in the layer. 

Korea horizontal tfornce in. “the x,-direction there is 
no displacement corresponding to m = 9. The displacement 
corresponding to m = 1 is: 

2 vtl1 
F J. (tr) Soe CE) a econ 
(SH) 2 il ) 1 2 -NBH 
= &s ———. K A.38 
(3) 27,8 IDE Ses er eneh d(Cr) eg] 2 b eR 
and if the receiver ison the x, ~axis, ¢ = 7/2 and the 


displacement is 
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Restoring the integration with respect to ¢ and 
performing the change of variable 
GC o= kx (A.40) 


we write the displacement of the S wave in the form 
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For-a horizontal force in the x,-direction and a 


(A.41) 


receiver on the x, axis, we get after similar steps: 


ee = ae. e where 
—(S) | = ; 1 ee =~ wea x tn, 
q 7 re ) | sree J, (k xrJe dx (A.42) 
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0 vx°4+n 
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b's Compressional-Shear Waves 


The system for the coefficients Ee Mpa wd te Can 
man? om mn. im am 
be brought to the following form 
A (H-d) eee e (na 
me ‘< ee SV 
Ba a4 ee e 3‘ ) (A.44) 


m x ~m 


where Av. is, a. (6X6). matrix, identical with that obtained 
~m 
taking ithe: first six columns,of theomatrix (23) in 


Abramovici (1970), 


Ei Bi /o,> A= k Vx +e (A.45) 


oh] being the P-velocity in the layer, and 
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2 +1 A(H-2d) 2 T 
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Weve Van 2 Cee 
7 x a Ve mes 
poo) = m+] ee B(H-2d) m+] eet B(H-2d) 
Bagi | (cl) — 3 o GI eo ieis e 2 (A.46) 
vx" +n ri 
i 
i 2 +e Bee BOs) | “ 
3 > y > 
a ae Veen? eae Sens 
aL 1 1 
Pemedning transposed, Thus, the displacement is split in a 
natural way into P and SV motions 
P SV 
oe eg u‘ 7 u‘ (A.47) 
m m m 
with Gay nd cee satisfying: 
m m 
) : = is 
UNS een (ERS poste Seon 118s) (A. 48) 
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Solving these systems of equations, expanding into 
generalized rays and taking into account the expression 
(A817). fory the tharmonic ivectons Aweuget (the (vertical, radial 
and angular displacement components wards and . es 
respectively for the P and SV motions as integrals with 


respect to x, having as integrands: 
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is 


SEP) 2 im ~AMH- 
ie -k, Vx te, J 6k xr) Re(i e %) ) R, (u,v) ee 
Usv 
=(P) 2 -~ im -AMH-— 
9 40k = -k)x Jk, xr) Re (ie >) ) S,(u,ve see 
Uyv 
7{P) =- "3 (k xr)Re(i ate ) S phase 
m Tem el m° : 3feeve 
u,v 
/ (A.49) 
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R (u,v), S$, (u,v) being defined in Abramovici (1970) and 
Ry (uy 5 S,(u,v) in Abramovici and Gal-Ezer (1978). 
For a horizontal force in the x,-direction, we have 


to take only m = 1. When the receiver is located on the 


x, ~axis, the displacement is: 
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For a horizontal force in the x.-direction we have 


1 
again only the term m = 1. When the receiver is on the 
x,~axis, the displacement is: 

Fe) @ (SV) 2 
F * xJ_(k, xr) 
ube gue a 21S ee oa Eg -AMH-B vH 
wW unk, ) u R,(u,ve dx 
LH» 0 
Z val SD ED 
F i te. vx +n 
-(SV)_ 1 1 3k R, GAUH-BNH (A.52) 
pw ss Tk, ) aprarrmerrn are epg sea NUT ROG 
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F x ath (k, xr) 
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A4. The solution in the time domain 


Using the generalized Cagniard-Pekeris method 
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(Abramovici, 1978), each term in (A.36) and (A.49) can be 
inverted directly into the time domain if the time varia- 
tion of the stress-discontinuity is a power of time t 
multiplying a step function or any linear combination of 
such functions. In particular, if we have a concentrated 


force Facting in’ the horizontaleplane (x »X,) of components 


1 


F(t) = £,U(t)t 
Caos) 


F(t) = f£,U(t)t 


where f,>f, ane scale factors and U(t) is the unit step 


function, we have 


eo ed eee (A.54) 


and the displacement components, according to GAS A 1) Chora) 


ands (CA'.51)-CA.52), L0F .2. receiver on tie x,-axis, are: 
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e f 3 
amp” ee 0 am 
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The inversion into the time domain is expressed in 
terms of some special quantities Ba ‘ (Abramovici, 1978) 


that depend on two functions and on the time variable: 
Fo? gee (bere (A.58) 
m m 


The function ¢ is the "amplitude" of the integrands in (A(55),(56), 
(57)) whereas g is the argument of the exponentials. For 
simplicity, we shall drop here the lower index. 
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horizontal force are: 
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p being the non-dimensional epicentral distance 
= 57H (A.64) 


and pop, the critical epicentral distance for a head wave to 
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occur. Besides, tT is a non-dimensional time parameter, 


To the arrival time of the geometrical ray reflected uz 
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corresponding head wave. The expressions for X ’ , Y ’ 
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For the particular cases needed in (A.59)," CA.60); 
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APPENDIX B 
In general, the steady-state plane wave solution to 


(3.3) is a constant times 


i eee _ AS 1 around 
ai Dat &) ae 3 Aer ee r-wt) (B41) 
where 
oe: Poa A (B22) 
> 
Peis the propagation vector” and A is the attenuation 
vector. 
From (B.2) we have 
Soles se Re[k* | (Bs) 
and 
Z 
2PAcosy = Im[K* ] (Bir 4) 


a 
where P and A are the amplitudes of P and A and Y- "Ls 
the angle between them. We now define the loss factor 


for P and S waves as follows: 


G = At2u for P waves 
nd ABLE cece (B.5) 
Re[G] 
fase =u for S waves. 
ay, i (B.6) 
Thus G = m(1+ iq) , 
and me = rie = Chi) 


where p is the density and m is a constant to be determined. 
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MBO), CB. 7) we get: 


P= w 2 VE+F (B.8) 
A = wo ye VE-F (B.9) 


2 
where F = SUT aaa CBee.) 
2(Q°+1) , 
and E = A + Re sree? ; (B.11) 
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But the phase velocity c(w) is equal to 


mn = e- (rn GEer) CB 2) 


Therefore using (B.6) the complex modulus G is obtained 


aS) LOL ows: 
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Futterman's model provides the phase velocity and Q as 


a Puoet1on of Trequency:. 


Equation (8.2) impiies that sin ihe, most peneran case 


p and A are not parallel. When y=0 the wave is called 


homogeneous while when y#0 the wave is called inhomogeneous. 


In our examples we assume that the source transmits 


homogeneous waves, i.e. y=0. Care must be taken working 
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out the angles of emergence for the reflecting waves. 
Consider a plane P wave striking an interface and producing 
reflected and transmitted P and S waves as shown in Figure 
2>.- For simplicity the transmitted P‘and S waves and the 
reflected S waves are not shown in Figure 25. If we know 
the angle of incidence oy and the attenuation angle vy 

then the angle of emergence 8. and the Y, can be calculated 
using sonell's Law, swhichystatesathatertheshorizontal wave 


number K. isi*constent Lockett. 1962,<Borenerdt. 197 /) 


i.e. Re[K ] = constant 
CBRL) 
or Pp, sing, = P,sin 6, 
Im[K,] = constant 
i - = i - : Biel 
or A,sin(6, Y,? A,sin(6, Y5) ( 5) 


Equations (B.14), (B.15) are used to calculate 6, and Yo: 
Similarly, one can calculate angles of transmission and 


emergence for the transmitted and reflected P and S waves. 
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List of Computer Programs 


THIS PROGRAM TAKES A GENERALIZED RAY AND 
INTRODUCES THE EFFECTS OF ANELASTICITY 


ODNOWAWNDN = 


(Seka ekaleh(eate (eee (ie eae rele) (a) lilo) @ite) @) le) @) 


100 


101 


ACCORDING TO FUTTERMAN’S MODEL 


V1=P-=VELOCITY 
V2=S=VEEOCLITY 
V3=P-VELOCITY 
V4=S-VELOCITY 
R1=DENSITY IN 
R2=DENSITY IN 


IN THE LAYER 

IN THE LAYER 

IN THE HALF SPACE 
IN THE HALF SPACE 
THE LAYER 

THE HALF SPACE 


Q1Z=Q VALUE IN THE LAYER 
Q2Z=Q VALUE IN THE HALF SPACE 
DEPTH=DEPTH-OF THE SOURCE 


DIST=HORIZONTAL DISTANCE SOURCE-RECEIVER 


DEL=SAMPLING RATE 
FL=LOW FREQUENCY 
FH=HIGH FREQUENCY 
GAM=ATTENUATION ANGLE 


X=GENERALIZED 


XR=GENERALIZED RAY WITH THE EFFECTS OF ANELASTICITY 
SUBROUTINE PEL CALCULATES REFLECTION COEFFICIENTS 


RAY 


FOR A P INCIDENT PLANE WAVE 


SUBROUTINE SEL CALCULATES REFLECTION COEFFICIENTS 


FOR AN S INCIDENT PLANE WAVE 


SUBROUTINE 
SUBROUTINE 
SUBROUTINE 
SUBROUTINE 
SUBROUTINE 


PVEL IS THE ANELASTIC VERSION OF PEL 
SVEL IS THE ANELASTIC VERSION OF SEL 
FASTF PERFORMS A FAST FOURIER TRANSFORM 
PHAMP PUTS THE PHASE BETWEEN O AND TWOPI 
LEQT1C SOLVES A SYSTEM OF FOUR EQUATIONS 


WITH FOUR UNKNOWNS. 


REAL X(1024) ,XR( 1024) ,X1( 1024) ,AAM( 1024) ,PPH( 1024) 


REAL AM( 1024), 
REAL VP1( 1024) ,VP2(1024),VS1( 1024) ,VS2( 1024) ,QP( 1024) ,Q5( 1024) 


REAL F(1024),Z(1024),17( 1024) ,FG(1024) 


PH( 1024) 


COMMON /AREA1/V1,V2,V3,V4,R1,R2 


COMMON /AREA2/VP1,VP2,VS1,VS2,QP,QS,RS1,RS2,GAM 


READ THE DATA,VELOCITIES,ELASTIC PARAMETERS ETC. 


READ(5,100)HLAY 
FORMAT(6X,F9.4) 
READ(5,100)DEPTH 
READ(5,100)DIST 


READ(5, 100)DEL 


READ(5,100)FL 
READ(5,100)FH 
READ(5,100)V1 
READ(5, 100) V2 
READ(5, 100) V3 
READ(5,100)V4 


READ(5,100)Q1Z 
READ(5, 100) Q2Z 


READ(5,100)R1 
READ(5, 100)R2 


READ(5,100)GAM 


READ(5, 101) INDEXP 


FORMAT(8X,12) 


READ(5,101)INDEXS 


READ(5,101)KK 
READ(5,101)KL 
READ(5,101)KM 
READ(5,101)KN 
RS1=R1 
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RS2=R2 
READ THE GENERALIZED RAY 

N=512 

N2=N/2 

PI=3.141593 

NX=1 

CONTINUE 

READ( 1,299, END=58)X(NX) 
FORMAT (40x ,D20.9) 

NX=NX+ 4 

GO TO 59 

CONTINUE 

NX=NX-1 

TAPERING AND ADDING ZEROS 
ND2=(N/2)+1 

NTA=476 

NEND=N-NTA 

NX 1=NX+4 

MTOT=N-NX-NTA 
RMTOT=FLOAT(MTOT) 

IC=0 

DO 2 I=NX1,NEND 

IC=IC+1 

X (I) =X (NX) *((1.0+COS(PI*FLOAT(IC)/RMTOT))*0.5) 
NENDS=NEND+ 1 

DO 3 I=NENDS,N 

X(1)=0.0 

FOURIER TRANSFORM OF THE RAY 
DO?7 (P=45N 

XR(1I)=X(1) 

XI(1I)=0.0 

CONTINUE 

CALL FASTF(XR,XI,N) 

DO 10 I=1,ND2 
AAM(I)=((XR( I) *XR( I) )+(X1 (1) *X1(1))) **0.5 
PPH(1)=ATAN2(XI(1I),XR(I)) 
KK=+1 DIRECT P OR S 

KK=-1 PP,PS,SP,SS 

KL=+1 DIRECT P 

KL=-1 DIRECT §S 

KM=O SP OR SS 

KM=+1 PP 

KM=-1 PS 

KN=+1 SS 

KN=-1 SP 

INDEXP=+1 PP 

INDEXP=+2 PS 

INDEXS=+1 SP 

INDEXS=+2 SS 

COMPUTE THE PATH TRAVELLED AND THE ANGLE OF INCIDENCE 
IF(KK.GT.O) GO TO 251 
RO=((((2.O*HLAY) -DEPTH) **2)+(DIST**2))**0.5 
VC=DIST/((2.0*HLAY ) -DEPTH) 
AN=ATAN( VC) 

RO2=HLAY/COS(AN) 

RO1=RO-RO2 

GO2.TO? 252 
RO=((DIST**2)+(DEPTH**2) )**O.5 
CONTINUE 
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510 


Silt 


332 


512 


331 


513 


334 


EC=0.5772 
XH=ALOG(FH/FL) 

VP 1L=V1*(1.0-(1.0/(P1*Q1Z) )*(EC+XH) ) 

VS 1L=V2*(1.0-(1.0/(P1I*Q2Z) )*(EC+XH) ) 
VP2L=V3*(1.0-(1.0/(P1*Q1Z) ) *(EC+XH) ) 
VS2L=V4*(1.0-(1.0/(P1*Q2Z) )*(EC+XH) ) 
VP1(1)=VP4L 

VS1(1)=VS4L 

VP2(1)=VP2L 

VS2(1)=VS2L 

QP(1)=Q12Z 

QS(1)=Q2Z 

DO 510 I=2,ND2 
F(I)=2.0*PI*FLOAT(I)/(FLOAT(N) *DEL) 
Z(1)=ALOG(F(1I)/FL) 
VP1(1)=VP1IL*(1.0-(1.0/(P1*Q1iZ) )*(EC+Z2(1)))**(-1) 
VS1(1)=VS1L*(1.0-(1.0/(P1*Q2Z) )*(EC+Z(1)))**(-1) 
VP2(I)=VP2L*(1.0-(1.0/(PI*Q1Z))*(EC+Z(1)))**(-1) 
VS2(1)=VS2L*(1.0-(1.0/(PI*Q2Z) )*(EC+Z(1I)))**(-1) 
QP(1)=Q12-((1.0/P1)*(EC+Z(1))) 
QS(1)=Q2Z-((1.0/P1)*(EC+Z(1))) 

CONTINUE 

INTRODUCE ATTENUATION AND DISPERSION ALONG THE PATH 
THE EFFECT OF INTERFACE IS ALSO CONSIDERED 
IF(KK.LT.O) GO TO 331 

FECKE LTO) “GOTO 332 

AM(1)=AAM(1) 

PH(1)=PPH(1) 

DO 511 I=2,ND2 
Di=F(1I)*RO*((1.0/V1)-(1.0/VP1(1))) 
PH(I)=PPH(1I)+D1 

AL=F(I)/(2.0*QP(1I)*VP1(I)) 

EK=-AL*RO 

AM(I)=AAM(1)*EXP(EK) 

GO TO 3900 

AM(1)=AAM( 1) 

PH(1)=PPH(1) 

DO 512 I=2,ND2 
D1=F(1)*RO*((1.0/V2)-(1.0/VS1(I))) 
PH(1I)=PPH(1I)+D1 

AL=F(1)/(2.0*QS(1)*VS1(I)) 

EK=-AL*RO 

AM(1I)=AAM(1) *EXP(EK) 

GO TO 900 

IF(KM.EQ.0) GO TO 333 

IF(KM.LT.O) GO TO 334 

CALL PEL(AN, INDEXP,AME,PHE) 

AM(1)=AAM(1) 

PH(1)=PPH(1) 

DO 513 I=2,ND2 
D1i=F(1I)*RO*((1.0/V1)-(1.0/VP1(1))) 

IM=I 

CALL PVEL(IM,AN, INDEXP, AMV, PHV) 
PH(1)=PPH(1I)+D1+PHV 
AL=F(1I)/(2.0*QP(1I)*VP1(1I)) 

EK=-AL*RO 

AM(1)=AAM(I) *EXP( EK) *AMV/AME 

GO TO 900 

AM(1)=AAM( 1) 

PH(1)=PPH(1) 
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CALL PEL(AN, INDEXP, AME, PHE) 

DO 514 I=2,ND2 
Di=FCI)*RO1* (C1. 0/V1)= (400/ VP 11 
D2=F(1)*RO2*((1.0/V2)-(1.0/VS1(I 
IM=I 

CALL PVEL(IM,AN, INDEXP, AMV, PHV) 
PH(1)=PPH(1)+D1+D2+PHV 
AL1=F(1I)/(2.0*QP(1)*VP1(I)) 
EK1=-AL1*RO14 
AL2=F(1I)/(2.0*QS(1)*VS1(T)) 
EK2=-AL2*RO2 

EK=EK1+EK2 

AM(1)=AAM(I1) *EXP(EK) *AMV/AME 

GO TO 3900 

IF(KN.LT.O) GO TO 335 
AM(1)=AAM(1) 

PH(1)=PPH(1) 

CALL SEL(AN, INDEXS,AME,PHE) 

DO 515 I=2,ND2 
D1=F(1)*RO*((1.0/V2)-(1.0/VS1(I1) 
IM=I 

CALL SVEL(IM,AN, INDEXS, AMV, PHV) 
PH(1)=PPH(1)+D1+PHV 
AL=F(1I)/(2.0*QS(1I)*VS1(I)) 
EK=-AL*RO 

AM(I)=AAM(1) *EXP(EK) *AMV/AME 

GO TO 800 

AM(1)=AAM(1) 

PH(1)=PPH(1) 

CALL SEL(AN, INDEXS, AME, PHE) 

DO 516 I=2,ND2 
Di=F(1)*RO1*((1.0/V2)-(1.0/VS1(I 
D2=F(1)*RO2*((1.0/V1)-(141.0/VP1(I 
IM=I 

CALL SVEL(IM,AN, INDEXS, AMV, PHV) 
PH(I)=PPH(1I)+D1+D2+PHV 
AL1=F(1)/(2.0*QS(1)*VS1(1)) 
EK1=-AL1*RO1 
AL2=F(1I)/(2.0*QP(1)*VP1(1)) 
EK2=-AL2*RO2 

EK=EK1+EK2 
AM(1I)=AAM(I) *EXP( EK) *AMV/AME 
CONTINUE 

INVERSE FOURIER TRANSFORM 

DO 681 I=1,ND2 
XR(1)=AM(1)*COS(PH(I)) 
XI(1)=AM(1)*SIN(PH(I)) 

CONTINUE 

DO 682 I=2,N2 
XR(ND2+1-1)=XR(ND2-I+1) 
XI(ND2+I-1)=-XI(ND2-I+1) 
CONTINUE 

DO 683 I=1,N 

XI(I)=-XI(1) 

CONTINUE 

XR(ND2)=0.0 

XI(ND2)=0.0 

CALL FASTF(XR,XI,N) 

DO 684 I=1,120 
XR(I)=XR(1I)/(FLOAT(N) ) 
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CONTINUE 

WRITE( 19,298) (XR(1),1=1,100) 
FORMAT (40X ,D20.9) 

STOP 

END 

SUBROUTINE FASTF(FR,FI,N) 


N IS THE NUMBER OF DATA POINTS = 2**M 

FR IS THE REAL DATA SET 

FI IS THE IMAGINARY PART OF THE DATA SET (= 0.0 IF ONLY REAL) 
FIRST COMPUTE M 


REAL FR(N), FI(N), GR, GI, ER, EI,EU, EZ 
M=0 

KD=N 

KD=KD/2 

M=M+ 1 

LECKO.. Ges 2) GO TO" 1 

ND2 = N/2 

NM1=N- 14 

L=1 


SHUFFLE INPUT DATA IN BINARY DIGIT REVERSE ORDER 


DO 4 K=1,NM1 
TRC GE a tb GO at 0-2 
GR=FR(L) 

GI=FI(L) 

FR(L)=FR(K) 

FI(L)=FI(K) 

FR(K)=GR 

FI(K)=GI 

NND2=ND2 

IF(NND2 .GE. L) GO TO 4 
L=L-NND2 

NND2=NND2/2 

GO TO 3 

L=L+NND2 

PI=3.14159265 


FIRST ARRANGE ACCOUNTING OF M STAGE 


DO 6 J=1,M 
NJ=2**J 
NUD2=NuU/2 
EU=1.0 

EZ=0.0 
ER=COS(-PI/NUD2) 
EI=SIN(-PI/NUD2) 


COMPUTE FOURIER TRANSFORM IN EACH M STAGE 


DO 6 IT=1,NJD2 

DO 5 IW=IT,N,NdJ 
IWJ=IW+NJD2 

GR=FR(IWJ) *EU-FI( IW) *EZ 
GI=FI(IWJ) *EU+FR(IWJ) *EZ 
FR(IWJ)=FR(IW)-GR 
FI(IWJ)=FI(IW)-GI 
FR(IW)=FR(IW)+GR 
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FIC IW)=FI(IW)+GI 

SEU=EU 

EU=SEU*ER-EZ2*EI 

EZ=EZ*ER+SEU*EI 

RETURN 

END 

SUBROUTINE PEL(D,INDEXP,AM,PH) 
REAL WA(4) 

COMPLEX HS,A(4,4),B(4,1),CABS,CSQRT 
COMMON /AREA1/V1,V2,V3,V4,R1,R2 
HS=SIN(D)/V1 

A(1,1)=-V1*HS 
A(1,2)=-CSQRT(1.0-((HS*V2)**2)) 
A(1,3)=V3*HS 
A(1,4)=CSQRT(1.0-((HS*V4) **2) ) 
A(2,1)=CSORT(1.0-((HS*V1)**2)) 
A(2,2)=-V2*HS 
A(2,3)=CSQRT(1.0-((HS*V3) **2)) 
A(2,4)=-V4*HS 
A(3,1)=2.0*R1*V2*V2*HS*A(2,1) 
A(3,2)=R1*V2*(1.0-2.0*V2*V2*HS*HS ) 
A(3,3)=2.0*R2*V4*V4*HS*A(2,3) 
A(3,4)=R2*V4*(1.0-2.0*V4*V4*HS*HS ) 
A(4,1)=-R1*V1*(1.0-2.0*V2*V2*HS*HS) 
A(4,2)=2.0*R1*V2*V2*HS*(-A(1,2)) 
A(4,3)=R2*V3*(1.0-2.0*V4*V4*HS*HS ) 
A(4,4)=-2.0*R2*V4*V4*HS*A(1,4) 
B.G1,.4) =a Ast o1,) 

B(2,1)=A(2,1) 

B(3,1)=A(3,1) 

B(4,1)=-A(4,1) 

CALL LEQT1C(A,4,4,B8,1,4,0,WA,IER) 
AM=CABS(B(INDEXP, 1)) 

CALL PHAMP(B(INDEXP,1),AM,PH) 
RETURN 

END 

SUBROUTINE PHAMP(POM, AMP, PHASE) 
COMPLEX POM,CABS 

AMP =CABS(POM) 

REZ=REAL(POM) 

AMZ=-AIMAG(POM) 

IF (AMP.EQ.0) GO TO 27 

IF (REZ.EQ.0) GO TO 28 

IF (AMZ.EQ.0) GO TO 29 
AR=ATAN2(AMZ,REZ) 

IF (AR.LT.O) GO TO 30 

PHASE=AR 

RETURN 

PHASE =AR+2*3. 141593 

GO TO 31 

PHASE=O 

GO TO 31 

IF (AMZ.GT.O) GO TO 32 

PHASE=3. 141593*1.5 

GO TO 31 

PHASE=3. 141593*0.5 

IF (REZ.GT.O) GO TO 27 

PHASE=3. 141593 

GO TO 31 

END 
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SUBROUTINE PVEL(NW,D, INDEXP,AM, PH) 


REAL WA(4),VP1( 1024) ,VvP2( 1024) ,VS1( 1024) ,VS2( 1024) 


REAL QP( 1024) ,QS( 1024) 
REAL MiR,M1I,M2R,M21,M3R,M3I,M4R,M4I 


COMPLEX HS,A(4,4),B(4,1),CABS,CSQRT,CMPLX 


COMPLEX M1,M2,M3,M4,VV1,VV2,VV3,VV4 


COMMON /AREA2/VP1,VP2,VS1,VS2,QP,QS,RS1,RS2,GAM 


PI=3.141593 
GA=P1*GAM/ 180.0 
$1=2.0*(1.0+(1.0/(QP(NW) **2))) 
S2=1.0+(1.0+(1.0/(QP(NW) **2)))**0.5 
M1R=(VP1(NW)**2)*RS1*S2/S1 

M11 =M1R/QP(NW) 

M1=CMPLX(M1R,M1I) 

VV 1=CSORT(M1/RS1) 
T1=2.0*(1.0+(1.0/(QS(NW)**2))) 
T2=1.0+(1.0+(1.0/(QS(NW) **2)))**O.5 
M2R=(VS1(NW)**2)*RS1*T2/T1 

M21 =M2R/QS(NW) 

M2=CMPLX(M2R,M21) 
VV2=CSOQRT(M2/RS1) 
TT1=2.0*(1.0+(1.0/(QP(NW)**2))) 
TT2=1.0+(1.0+(1.0/(QP(NW) **2)))**0.5 
M3R=(VP2(NW)**2)*RS2*TT2/TT1 
M3I=M3R/QP (NW) 

M3=CMPLX(M3R,M31) 
VV3=CSQRT(M3/RS2) 
SS1=2.0*(1.0+(1.0/(QS(NW) **2))) 
$S2=1.0+(1.0+(1.0/(QS(NW) **2)) )**0.5 
M4R=(VS2(NW) **2)*RS2*SS2/SS1 

M41 =M4R/QS(NW) 

M4=CMPLX(M4R,M41) 
VV4=CSQRT(M4/RS2) 
X1=1.0+(1.0/(QP(NW) *COS(GA)))**2 
X11=SQRT(X1) 
X2=1.0+(1.0/(QP(NW)))**2 
X22=SQRT(X2) 

XN=1.0+X14 

XD=1.0+X22 

XZ=XN/XD 

XV=SORT(XZ) 

PP1=XV/VP1(NW) 

XNN=X11-1.0 

XZZ=XNN/XD 

XVV=SQRT(XZZ) 

AA1=XVV/VP1(NW) 

DD1=D-GA 

HSR=SIN(D)*PP1 

HSI=-SIN(DD1)*AA1 
HS=CMPLX(HSR,HSI) 

A(1,1)=-VV1*HS 
A(1,2)=-CSQRT(1.0-((HS*VV2)**2)) 
A(1,3)=VV3*HS 

A( 1,4) =CSQRT(1.0-((HS*VV4)**2)) 
A(2,1)=CSQRT(1.0-((HS*VV1)**2)) 
A(2,2)=-VV2*HS 
A(2,3)=CSQRT(1.0-((HS*VV3)**2)) 
A(2,4)=-VV4*HS 

A(3, 1)=2.0*RS1*VV2*VV2*HS*A(2, 1) 
A(3,2)=RS1*VV2*(1.0-2.0*VV2*VV2*HS*HS ) 
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A(3,3)=2.0*RS2*VV4*VV4*HS*A(2,3) 
A(3,4)=RS2*VV4*(1.0-2.0*VV4*VV4*HS*HS ) 
A(4,1)=-RS1*VV1*(1.0-2.0*VV2*VV2*HS*HS) 
A(4,2)=2.0*RS1*VV2*VV2*HS*(-A(1,2)) 
A(4,3)=RS2*VV3*(1.0-2.0*VV4*VV4*HS*HS ) 
A(4,4)=-2.0*RS2*VV4*VV4*HS*A(1,4) 
B(1,1)=-A(1,1) 

B( 2.0%). =A (2, 1) 

B(3.<1),=A.63.,.1.) 

B(4,1)=-A(4,1) 

CALL LEQT1C(A,4,4,B,1,4,0,WA,IER) 
AM=CABS(B(INDEXP, 1)) 

CALL PHAMP(B(INDEXP,1),AM,PH) 

RETURN 

END 

SUBROUTINE SEL(D, INDEXS,AM, PH) 

REAL WA(4) 

COMPLEX HS,A(4,4),B(4,1),CABS,CSQRT 
COMMON /AREA1/V1,V2,V3,V4,R1,R2 
HS=SIN(D)/V1 

A(1,1)=-V1*HS 
A(1,2)=-CSQRT(1.0-((HS*V2) **2) ) 
A(1,3)=V3*HS 
A(1,4)=CSQRT(1.0-((HS*V4) **2)) 
A(2,1)=CSQRT(1.0-((HS*V1)**2)) 
A(2,2)=-V2*HS 
A(2,3)=CSQRT(1.0-((HS*V3) **2) ) 
A(2,4)=-V4*HS 
A(3,1)=2.0*R1*V2*V2*HS*A(2, 1) 
A(3,2)=R1*V2*(1.0-2.0*V2*V2*HS*HS ) 
A(3,3)=2.0*R2*V4*V4*HS*A(2,3) 
A(3,4)=R2*V4*(1.0-2.0*V4*V4*HS*HS ) 
A(4,1)=-R1*V1*(1.0-2.0*V2*V2*HS*HS) 
A(4,2)=2.0*R1*V2*V2*HS*(-A(1,2)) 
A(4,3)=R2*V3*(1.0-2.0*V4*V4*HS*HS ) 
A(4,4)=-2.0*R2*V4*V4*HS*A(1,4) 

B(1, 1}#-A(1,.2) 

B(2,1)=A(2,2) 

B(3,1)=A(3,2) 

B(4,1)=-A(4,2) 

CALs LEQTACGCA4. 4,.B.154.0,WA,LER) 
AM=CABS(B(INDEXS, 1) ) 

CALL PHAMP(B(INDEXS,1),AM,PH) 

RETURN 

END 

SUBROUTINE SVEL(NW,D, INDEXS,AM,PH) 
REAL WA(4),VP1(1024),VP2( 1024) ,VS1(1024),VS2( 1024) 
REAL QP(1024),QS( 1024) 

REAL MiR,M11,M2R,M2I] ,M3R,M3I,M4R,M4I 
COMPLEX HS,A(4,4),B(4,1),CABS,CSQRT ,CMPLX 
COMPLEX M1,M2,M3,M4,VV1,VV2,VV3, VV4 
COMMON /AREA2/VP1,VP2,VS1,VS2,QP,QS,RS1,RS2,GAM 
PI=3.141593 

GA=PI*GAM/ 180.0 
S1=2.0*(1.0+(1.0/(QP(NW)**2))) 
S2=1.0+(1.0+(1.0/(QP(NW) **2)))**0.5 
MiR=(VP1(NW) **2)*RS1*S2/S1 
MiI=M1R/QP(NW) 

M1=CMPLX(MiR,M11) 

VV1=CSQRT(M1/RS1) 
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T1=2.0*(1.0+(1.0/(OS(NW) **2))) 
T2=1.0+(1.0+(1.0/(OS(NW) **2)))**0.5 
M2R=(VS1(NW) **2)*RS1*T2/T1 
M21=M2R/QS(NW) 

M2=CMPLX(M2R,M21) 

VV2=CSQRT(M2/RS1) 

TT 1=2.0*(1.0+(1.0/(OP(NW) **2))) 
TT2=1.04+(1.0+(1.0/(QP(NW) **2)))**0.5 
M3R=(VP2(NW)**2)*RS2*TT2/TT1 

M31 =M3R/QP(NW) 

M3=CMPLX(M3R,M31) 

VV3=CSORT(M3/RS2) 
$S1=2.0*(1.0+(1.0/(QS(NW) **2))) 
S$S2=1.0+(1.0+(1.0/(QS(NW) **2)))**0.5 
M4R=(VS2(NW) **2) *RS2*SS2/SS1 
M41=M4R/QS(NW) 

M4=CMPLX(M4R,M41) 

VV4=CSORT(M4/RS2) 
X1=1.0+(1.0/(QP(NW) *COS(GA)))**2 
X11=SORT(X1) 
X2=1.0+(1.0/(QP(NW)))**2 
X22=SQRT(X2) 

XN=1.0+X11 

XD=1.0+X22 

XZ=XN/XD 

XV=SOQRT(XZ) 

PP1=XV/VP1(NW) 

XNN=X11-1.0 

XZZ=XNN/XD 

XVV=SQRT(XZZ) 

AA1=XVV/VP1(NW) 

DD1=D-GA 

HSR=SIN(D)*PP 1 

HSI=-SIN(DD1)*AA1 

HS=CMPLX(HSR,HSI) 

A(1,1)=-VV1*HS 
A(1,2)=-CSQRT(1.0-((HS*VV2) **2)) 
A(1,3)=VV3*HS 
A(1,4)=CSQRT(1.0-((HS*VV4) **2) ) 
A(2,1)=CSQRT(1.0-((HS*VV1)**2) ) 
A(2,2)=-VV2*HS 
A(2,3)=CSQRT(1.0-((HS*VV3) **2) ) 
A(2,4)=-VV4*HS 
A(3,1)=2.0*RS1*VV2*VV2*HS*A(2, 1) 
A(3,2)=RS1*VV2*(1.0-2.0*VV2*VV2*HS*HS ) 
A(3,3)=2.0*RS2*VV4*VV4*HS*A(2,3) 
A(3,4)=RS2*VV4*(1.0-2.0*VV4*VV4*HS*HS ) 
A(4,1)=-RS1*VV1*(1.0-2.0*VV2*VV2*HS*HS ) 
A(4,2)=2.0*RS1*VV2*VV2*HS*(-A(1,2)) 
A(4,3)=RS2*VV3*(1.0-2.0*VV4*VV4*HS*HS ) 
A(4,4)=-2.0*RS2*VV4*VV4*HS*A(1,4) 
B(1,1)=-A(1,2) 

B(2,1)=A(2,2) 

B(3,1)=A(3,2) 

B(4,1)=-A(4,2) 

CALL LEQT1C(A,4,4,B,1,4,0,WA,IER) 
AM=CABS(B(INDEXS,1)) 

CALL PHAMP(B(INDEXS,1),AM,PH) 

RETURN 

END 
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PHASE UNWRAPPING PROGRAM 
THIS PROGRAM IS A MODIFIED VERSION OF 
TRIBOLET’S ORIGINAL PROGRAM 
REFERENCE: TRIBOLET, JU.M.( 1977). 
A NEW PHASE UNWRAPPING ALGORITHM. 
PEEE ;VOL) pASPP-25, NO 2 P. 170-177 
X=THE INPUT SIGNAL 

INITIALIZATION 


be2*r512 
M=9 
N=2**M 
N2=N/2 
ND2=N2+1 
PI=3.141596 
TWOPI=PI*2.0 
PARAMETERS FOR UNWRAPPED PHASE 


H=TWOPI/1024 
H1=H/L 
THLD1=1.8*PI 
THLD2=1.0*PI 
ISNX=+1 


TRANSFORM X(N) 
LOAD CXE AND CXO ARRAYS. 


DO 7 I=1,N 
CXE(1)=X(1) 
DO 8 I=1,N 
CxO(1I)=0.0 
CALL FASTF(CXE,CXO,N) 


TRANSFORM N*X(N) 
LOAD YR AND YI ARRAYS 


DO 19 I=1,N 
YR(I)=1*X(1) 

DO 20 I=1,N 
YI(1)=0.0 

CALL FASTF(YR,YI,N) 


CHECK IF SIGN REVERSAL IS REQUIRED 
IF(CXE(1).LT.0O) ISNX=-1 


COMPUTE LOGMAGNITUDE:STORE IN CXE 

COMPUTE PHASE DERIVATIVE:STORE IN YR 

COMPUTE PHASE PRINCIPAL VALUE:STORE IN YI 
FOR W=(TWOPI/N)*I,1=0,1,...N2 


DVTMN=O. 
DO 21 I=1,ND2 
A=CXE(I) 
B=CX0(1) 
C=YR(1) 

D =YI(I) 
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F=A*A+B*B 
YR(1)=-((A*C+B*D)/E) 
DVTMN=DVTMN+YR(T) 


CXE(1)=DLOG(E)/2 
IF(ISNX.EQ.+1)YI(I)=DATAN2(B,A) 
IF(ISNX.EQ.-1)YI(1I)=DATAN2(-B,-A) 
CONTINUE 
DVTMN=(2*DVTMN-YR(1)-YR(ND2))/N 
ELPINC=2*DVTMN*H 


PHASE UNWRAPPING BY ADAPTIVE INTEGRATION 


PH=O. 
DO 22 I=2,ND2 


STORE UNWRAPPED PHASE IN CXO 


CxO(I-1)=PH 


FORM PHASE ESTIMATE AT W=(TWOPI/N)*(I-1) BY 
TRAPEZOIDAL INTEGRATION WITH STEP SIZE TWOPI/N 


PHAINC=H*(YR(1)+YR(I-1)) 
IF (ABS(PHAINC-ELPINC).GT.THLD1) GO TO 23 
PH=CXO(1I-1)+PHAINC 


CHECK CONSISTENCY OF ESTIMATE 


AO=(PH-YI(1))/TWOPI 
A1=IFIX(AO)*TWOPI+YI(I) 
A2=A1+SIGN(TWOPI,AO) 

A3=ABS(A1-PH) 

A4=ABS(A2-PH) 
IF(A3.GT.THLD2.AND.A4.GT.THLD2) GO TO 23 


PHASE ESTIMATE WAS CONSISTENT 


PH=A1 

IF(A3.GT.A4)PH=A2 

IF (ABS(PH-CXO(I-1)).GT.P1) GO TO 23 
GO TO 22 


PHASE ESTIMATE WAS NOT CONSISTENT: 
ADAPT eSTERSSTZE 


INITIATE SOFTWARE STACK 


SP=1 
ISK(1)=N+1 

SK1(1)=YI(1) 
SK2(1)=YR(I) 


INITIATE REGISTERS 


IB=1 
B1i=CXO(I-1) 
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B2=YR(I-1) 


IF SOFTWARE STACK DIMENSION DOES NOT 
ALLOW FURTHER STEP REDUCTION, STOP 


24 IF((ISK(SP)-IB).GT.1) GO TO 25 
STOP 999 


DEFINE INTERMEDIATE FREQUENCIES (1.F) 
W=(TWOPI/N) *(1-2+(K-1)/N) 


25 K=(ISK(SP)+IB)/2 
COMPUTE DFTS OF X(N) AND N*X(N) AT I.F. 


AO=TWOPI*(I-2.+(K-1)/FLOAT(L))/FLOAT(N) 
CAH=( OCHO.) 
C2=(0. ,O. ) 
DO 26 J=1,N 
ARG=AO*(J-1) 
CO=CMPLX(COS(ARG) ,-SIN(ARG) ) *X(u) 
C1=C1+CO 
26 C2=C2+CO*(JU-1) 


COMPUTE PHASE DERIVATIVE AND THE 
PRINCIPAL VALUE OF THE PHASE AT I1.F. 


SP=SP+1 
ISK(SP)=K 

A=REAL(C1) 

B=AIMAG(C1) 

C=REAL(C2) 

D=AIMAG(C2) 

IFCISNX.EQ.+1) SK1(SP)=DATAN2(B,A) 
IF(ISNX.EQ.-1) SK1(SP)=DATAN2(-B,-A) 
SK2(SP)=-(A*C+B*D)/(A*A+B*B) 


EVALUATE ESTIMATE AT I.F. 


27 DELTA=H*(ISK(SP)-IB) 
PHAINC=DELTA*(B2+SK2(SP) ) 
IF (ABS(PHAINC-DELTA*2*DVTMN) .GT.THLD1) GO TO 24 
PH=B 1+PHAINC 


CHECK CONSISTENCY OF ESTIMATE AT I.F. 


AO=(PH-SK1(SP))/TWOPI 
A1=IFIX( AO) *TWOPI+SK1(SP) 
A2=A1+SIGN(TWOPI,AO) 

A3=ABS(A1-PH) 

A4=ABS(A2-PH) 
IF(A3.LT.THLD2.0R.A4.LT.THLD2) GO TO 28 


ESTIMATE WAS NOT CONSISTENT:REDUCE STEP SIZE 


GO TO 24 
ESTIMATE WAS CONSISTENT: UPDATE REGISTERS 


28 PH=A1 
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IF(A3.GT.A4)PH=A2 

IF (ABS(PH-B1).GT.PI) GO TO 24 
IB=ISK(SP) 

Bi=PH 

B2=SK2(SP) 

SP=SP-1 


WHEN SOFTWARE STACK IS EMPTY ,THE UNWRAPPED 
PHASE AT W=TWOPI*(I-1)/N IS HELD IN THE 
Bi REGISTER 


IF(SP.NE.O) GO TO 27 
PH=B1 


END OF STEP SIZE ADAPTATION 


CONTINUE 
IF(ISNX.EQ.-1)GO TO 191 
GO TO 192 

CONTINUE 

DO 193 I=1,N2 
CxO(1)=CXO(1)-PI 
CONTINUE 

STOP 

END 
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